.comp-ph] 24 Oct 2010 


ics 


1010.4992v1 [phys 


. 
. 


arXiv 


Applications of quantum Monte Carlo methods in condensed systems 


Jindřich Kolorenč 


Institute of Physics, Academy of Sciences of the Czech Republic, Na Slovance 2, 


18221 Praha 8, Czech Republic 


I. Institut für Theoretische Physik, Universität Hamburg, Jungiusstraße 9, 


20355 Hamburg, Germany 


kolorenc@fzu.cz 


Lubos Mitas 


Department of Physics and Center for High Performance Simulation, North Carolina 


State University, Raleigh, North Carolina 27695, USA 


lmitas@unity.ncsu.edu 


The quantum Monte Carlo methods represent a powerful and broadly applicable computa- 
tional tool for finding very accurate solutions of the stationary Schrodinger equation for 
atoms, molecules, solids and a variety of model systems. The algorithms are intrinsically 
parallel and are able to take full advantage of the present-day high-performance computing 
systems. This review article concentrates on the fixed-node/fixed-phase diffusion Monte 
Carlo method with emphasis on its applications to electronic structure of solids and other 


extended many-particle systems. 


PACS numbers: 02.70.Ss, 71.15.—m, 31.15.—p 


To appear in Rep. Prog. Phys. 


1 Introduction 
1.1 Many-body stationary Schrödinger equation . . . 


2 Methods a 
2.1 Variational Monte Carlo . . . o oaoa aaa 
2.2 Diffusion Monte Carlo ............... 

2.2.1 Fixed-node/fixed-phase approximation . . [E] 
2.2.2 Sampling the probability distribution El 
2.2.3 General expectation values ........ El 
2.2.4 Spin degrees of freedom .......... p 
2.3 Pseudopotentials . . .. o.oo pl 


3 From a finite supercell to the thermodynamic 


limit 
3.1 Twist-averaged boundary conditions ....... 
3.2 Ewald formula .. aoaaa 
3.3 Extrapolation to the thermodynamic limit . . .. 
3.4 An alternative model for Coulomb interaction 
ENErEY ea kee Be ee be ee E E e ed 


4 Trial wave functions 
4.1 Elementary properties ............... 
4:2: -Jastrow factor: os s m ac. e wa oa ce a ee Ga 
4.3 Slater—Jastrow wave function ........... 
4.4 Antisymmetric forms with pair correlations 
4.5 Backflow coordinates................ 


5 Applications 
5.1 Properties of the homogeneous electron gas 
5.2 Cohesive energies of solids ............. 
5.3 Equations of state .............200. 
5.4 Phase transitions ...............0.0., 
5.5 Lattice defects 
5.6 Surface phenomena................. 


5.7 Excited states... 2... 0.2... 0020084 
5.8 BCS-BEC crossover .............-00.% 
6 Concluding remarks 


1 Introduction 


Many properties of condensed matter systems can be 
calculated from solutions of the stationary Schrödinger 
equation describing interacting ions and electrons. The 
grand challenge of solving the Schrödinger equation 
has been around from the dawn of quantum mechanics 
and remains at the forefront of the condensed matter 
physics today and, undoubtedly, for many decades to 
come. 


The task of solving the Schrödinger equation for 
systems of electrons and ions, and predicting the quan- 
tities of interest such as cohesion and binding energies, 
electronic gaps, crystal structures, variety of magnetic 
phases or formation of quantum condensates is noth- 
ing short of formidable. Paul Dirac recognized this 
state of affairs already in 1929: “The general theory of 
quantum mechanics is now almost complete ... The 
underlying physical laws necessary for the mathemat- 


ical theory of a large part of physics and the whole 
chemistry are thus completely known, and the difficulty 
is only that the exact application of these laws leads 
to equations much too complicated to be soluble.” 
Arguably, this is the most fundamental approach to 
the physics of condensed matter: Applications of the 
rigorous quantum laws to models that are as close to 
reality as currently feasible. 

The goal of finding accurate solutions for stationary 
quantum states is hampered by a number of difficulties 
inherent to many-body quantum systems: 


(i) Even moderately sized model systems contain 
anywhere between tens to thousands of quantum 
particles. Moreover, we are often interested in 
expectation values in the thermodynamic limit 
that is usually reached by extrapolations from 
finite sizes. Such procedures typically require 
detailed information about the scaling of the 
quantities of interest with the system size. 


(ii) Quantum particles interact and the interactions 
affect the nature of quantum states. In many 
cases, the influence is profound. 


(iii) The solutions have to conform to quantum sym- 
metries such as the fermionic antisymmetry 
linked to the Pauli exclusion principle. This is 
a fundamental departure from classical systems 
and poses different challenges which call for new 
analytical ideas and computational strategies. 


(iv) For meaningful comparisons with experiments, 
the required accuracy is exceedingly high, espe- 
cially when comparing with precise data from 
spectroscopic and low-temperature studies. 


In the past, the most successful approaches to ad- 
dress these challenges were based mostly on reduction- 
ist ideas. The problem is divided into the dominant 
effects, which are treated explicitly, and the rest, which 
is then dealt with by approximate methods based on 
variety of analytical tools: perturbation expansions, 
mean-field methods, approximate transformations to 
known solutions, and so on. The reductionist ap- 
proaches have been gradually developed into a high 
level of sophistication and despite their limitations, 
they are still the most commonly used strategies in 
many-body physics. 

The progress in computer technology has opened a 
new avenue for studies of quantum (and many other) 
problems and has enabled researchers to obtain re- 
sults beyond the scope of analytic many-body theories. 
The performance of current large computers makes 
computational investigations of many-body quantum 


systems viable, allowing predictions that would be dif- 
ficult or impossible to make otherwise. The quantum 
Monte Carlo (QMC) methods described in this review 
provide an interesting illustration of what is currently 
possible and how much the computational methods 
can enrich and make more precise our understanding 
of the quantum world. 


Some of the ideas used in QMC methods go back 
to the times before the invention of electronic comput- 
ers. Already in 1930s Enrico Fermi noticed similarities 
between the imaginary time Schrodinger equation and 
the laws governing stochastic processes in statistical 
mechanics. In addition, based on memories of his col- 
laborator Emilio Segre, Fermi also envisioned stochas- 
tic methodologies for solving the Schrodinger equa- 
tion, which were very similar to concepts developed 
decades later. These Fermi’s ideas were acknowledged 
by Metropolis and Ulam in a paper from 1949 [2], 
where they outlined a stochastic approach to solv- 
ing various physical problems and discussed merits of 
“modern” computers for its implementation. In fact, 
this group of scientists at the Los Alamos National Lab- 
oratory attempted to calculate the hydrogen molecule 
by a simple version of QMC in the early 1950s, around 
the same time when a pioneering work on the first 
Monte Carlo study of classical systems was published 
by Metropolis and coworkers [3]. In the late 1950s, 
Kalos initiated development of QMC simulations and 
methodologies for few-particle systems and laid down 
the statistical and mathematical foundations of the 
Green’s function Monte Carlo method [4]. Eventually, 
simulations of large many-particle systems became 
practicable as well. First came studies of bosonic 
fluids modelling “He [5H7], and later followed investi- 
gations of extended fermionic systems exemplified by 
liquid *He [8] [9] and by the homogeneous electron gas 
{10} [71]. Besides these applications to condensed mat- 
ter, essentially the same methods were in mid-seventies 
introduced in quantum chemistry to study small molec- 
ular systems [12] [13]. To date, various QMC methods 
were developed and applied to the electronic structure 
of atoms, molecules and solids, to quantum lattice 
models, as well as to nuclear and other systems with 
contributions from many scientists. 


The term “quantum Monte Carlo” covers several 
related stochastic methodologies adapted to determine 
ground-state, excited-state or finite-temperature equi- 
librium properties of a variety of quantum systems. 
The word “quantum” is important since QMC ap- 
proaches differ significantly from Monte Carlo methods 
for classical systems. For an overview of the latter see 
for instance [14]. QMC is not only a computational 
tool for large-scale problems, but it also encompasses a 


substantial amount of analytical work needed to make 
such calculations feasible. QMC simulations often uti- 
lize results of the more traditional electronic structure 
methods in order to increase efficiency of the calcu- 
lations. These ingredients are combined to optimally 
balance the computational cost with achieved accuracy. 
The key point for gaining new insights is an appro- 
priate analysis of the quantum states and associated 
many-body effects. It is typically approached itera- 
tively: Simulations indicate the gaps in understanding 
of the physics, closing these gaps is subsequently at- 
tempted and the improvements are assessed in the next 
round. Such a process involves construction of zero- 
or first-order approximations for the desired quantum 
states, incorporation of new analytical insights, and 
development of new numerical algorithms. 


QMC methods inherently incorporate several types 
of internal checks, and many of the algorithms used 
possess various rigorous bounds, such as the varia- 
tional property of the total energy. Nevertheless, the 
coding and numerical aspects of the simulations are 
not entirely error-proof and the obtained results should 
be verified independently. Indeed, it is a part of the 
modern computational-science practice that several 
groups revisit the same problem with independent 
software packages and confirm or challenge the results. 
“Biodiversity” of the available QMC codes on the scien- 
tific market (including QWalk [15], QMCPACK [16], 
CHAMP [I7], CASINO [18], QMcBeaver and oth- 
ers) provides the important alternatives to verify the 
algorithms and their implementations. This is clearly a 
rather labourious, slow and tedious process, neverthe- 
less, experience shows that independently calculated 
results and predictions eventually reach a consensus 
and such verified data become widely used standards. 


In this overview we present QMC methods that 
solve the stationary Schrodinger equation for con- 
densed systems of interacting fermions in continuous 
space. Conceptually very straightforward is the varia- 
tional Monte Carlo (VMC) method, which builds on 
explicit construction of trial (variational) wave func- 
tions using stochastic integration and parameter opti- 
mization techniques. More advanced approaches rep- 
resented by the diffusion Monte Carlo (DMC) method 
are based on projection operators that find the ground 
state within a given symmetry class. Practical versions 
of the DMC method for a large number of particles 
require dealing with the well-known fermion sign prob- 
lem originating in the antisymmetry of the fermionic 
wave functions. The most commonly used approach to 
overcome this fundamental obstacle is the fixed-node 
approximation. This approximation introduces the 
so-called fixed-node error, which appears to be the key 


limiting factor in further increase in accuracy. As we 
will see in section [5] the fixed-node error is typically 
rather small and does not hinder calculation of robust 
quantities such as cohesion, electronic gaps, optical ex- 
citations, defect energies or potential barriers between 
structural conformations. By robust we mean quanti- 
ties which are of the order of tenths of an electronvolt 
to several electronvolts. Nevertheless, the fixed-node 
errors can bias results for more subtle phenomena, 
such as magnetic ordering or effects related to super- 
conductivity. Development of strategies to alleviate 
such biases is an active area of research. 

Fixed-node DMC simulations are computationally 
rather demanding when compared to the mainstream 
electronic structure methods that rely on mean-field 
treatment of electron-electron interactions. On the 
other hand, QMC calculations can provide unique in- 
sights into the nature of quantum phenomena and 
can verify many theoretical ideas. As such, they can 
produce not only accurate numbers but also new un- 
derstanding. Indeed, QMC methodology is very much 
an example of “it from bit” paradigm, alongside, for 
example, the substantial computational efforts in quan- 
tum chromodynamics, which not only predict hadron 
masses but also contribute to the validation of the 
fundamental theory [20] BI]. Just a few decades ago it 
was difficult to imagine that one would be able to solve 
the Schrödinger equation for hundreds of electrons by 
means of an explicit construction of the many-body 
wave function. Today, such calculations are feasible 
using available computational resources. At the same 
time, there remains more to be done to make the 
methods more insightful, more efficient, and their ap- 
plication less labourious. We hope that this review 
will contribute to the growing interest in this rapidly 
developing field of research. 

The review is organized as follows: The remainder 
of this section provides mostly definitions and nota- 
tions. Section [2] follows with description of the VMC 
and DMC methods. The strategies for calculation of 
quantities in the thermodynamic limit are presented 
in section [3] Section[4]introduces currently used forms 
of the trial wave functions and their recently devel- 
oped generalizations. The overview of applications 
presented in section [5] is focused on QMC calculations 
of a variety of solids and related topics. 


1.1 Many-body stationary Schrödinger 
equation 


Let us consider a system of quantum particles such as 
electrons and ions interacting via Coulomb potentials. 
Since the masses of nuclei and electrons differ by three 
orders of magnitude or more, the problem can be sim- 


plified with the aid of the Born—Oppenheimer approxi- 
mation, which separates the electronic degrees of free- 
dom from the slowly moving system of ions. The elec- 
tronic part of the non-relativistic Born—Oppenheimer 
hamiltonian is given by 


x 1 ZI A 1 
BS Ns Sine pa art (1) 


where i and j label the electrons and J runs over the 
ions with charges Zr. Throughout the review we em- 
ploy the atomic units, Me = h = e = 4r£o = 1, where 
Me is the electron mass, —e is the electron charge 
and £ọ is the permittivity of a vacuum. We are inter- 
ested in eigenstates |W,,) of the stationary Schrödinger 
equation 


A\W,) = En|Vn). (2) 


Colloquially, we call such solutions (either exact or 
approximate) and derived properties collectively the 
electronic structure. 

An important step forward in calculations of the 
eigenstates was made by Hartree [22] and Fock [23] by 
establishing the simplest antisymmetric wave functions 
and by formulating the Hartree-Fock (HF) theory, 
which correctly takes into account the Pauli exclu- 
sion principle [25]. The HF theory replaces the 
hard problem of many interacting electrons with a 
system of independent particles in an effective, self- 
consistent field. The theory was further developed by 
Slater [26] and others, and it has become a starting 
point of many sophisticated approaches to fermionic 
many-body problems. 

For periodic systems, the effective free-electron the- 
ory and the band theory of Bloch [27] were the first 
crucial steps towards our present understanding of 
the real crystals. In 1930s, Wigner and Seitz 
performed the first quantitative calculations of the elec- 
tronic states in sodium metal. Building upon the homo- 
geneous electron gas model, the density-functional the- 
ory (DFT) was invented by Hohenberg and Kohn [30] 
and further developed by Kohn and Sham who 
formulated the local density approximation (LDA) for 
the exchange-correlation functional. These ideas were 
later elaborated by including spin polarization [32], by 
constructing the generalized gradient approximation 
(GGA) [33] [34], and by designing a variety of orbital- 
dependent exchange-correlation functionals 37). 
The DFT has proved to be very successful and has 
become the mainstream computational method for 
many applications, which cover not only solids but 
also molecules and even nuclear and other systems 
[39]. The density-functional theory together with 
the Hartree-Fock and post-Hartree—Fock methods 


are relevant for our discussion of quantum Monte Carlo 
methodology, since the latter uses the results of these 
approaches as a reference and also for construction 
of the many-body wave functions. Familiarity with 
the basic concepts of the Hartree-Fock and density- 
functional theories is likely to make the subsequent 
sections easier to follow, but we believe that it is not 
a necessary prerequisite for understanding our exposi- 
tion of the QMC methods and their foundations. 


2 Methods 
2.1 Variational Monte Carlo 


In the variational Monte Carlo method, the ground 
state of a hamiltonian H is approximated by some trial 
wave function |W), whose form is chosen following a 
prior analysis of the physical system being investigated. 
Functional forms relevant to solid-state applications 
will be discussed later in section[4] Typically a number 
of parameters are introduced into |W), and these pa- 
rameters are varied to minimize the expectation value 
Egz = (Vr|H|Wy)/(Wr|Vr) in order to bring the 
trial wave function as close as possible to the actual 
ground state |W). 

Wave functions of interacting systems are non- 
separable, and the integration needed to evaluate Ey2 
is therefore a difficult task. Although it is possible to 
write these wave functions as linear combinations of 
separable terms, this tactic is viable only for a limited 
number of particles, since the length of such expansions 
grows very quickly as the system size increases. The 
variational Monte Carlo method employs a stochastic 
integration that can treat the non-separable wave func- 
tions directly. The expectation value Ey2, is written 
as 


|Wr(R)|? [AYr](R) 


pea N 
Ea =) Torin) mR E P 
N A 
~ Evmc 7 ee (3) 


where R = (11,12,...,7N) is a 3N-dimensional vec- 
tor encompassing the coordinates of all N particles 
in the system and the sum runs over M such vec- 
tors {R;} sampled from the multivariate probability 
density p(R) = |Wr(R)|?/(Wr|Ur). The summand 
EL (R) = [HV 7] (R)/Ur(R) is usually referred to as 
the local energy. We assume spin-independent hamil- 
tonians, and therefore spin variables do not explic- 
itly enter the evaluation of the expectation value 
This statement is further corroborated in section 
where the elementary properties of the trial wave func- 
tions |W) are discussed. 


Equation transforms the multidimensional in- 
tegration into a problem of sampling a complicated 
probability distribution. The samples {R;} can be 
obtained such that they constitute a Markov chain 
with transitions Ri+ı <— Ri; governed by a stochastic 
matrix M (Ri+ı +} Ri) whose stationary distribution 
coincides with the desired probability density p(R), 


pR’) = J M(R' 4 R)p(R) BYR forall R’. (4) 


After a period of equilibration, the members of the 
Markov sequence sample the stationary distribution 
regardless of the starting point of the chain, provided 
the matrix M(R’ 4+ R) is ergodic. Inspired by the 
way the samples explore the configuration space, they 
are often referred to as walkers. 

The Markov chain can be conveniently constructed 
with the aid of the Metropolis method [3] ÆI]. 
The transition matrix is factorized into two parts, 
M(R! + Ri) = T(R! + RiJ A(R’ + Ri), which cor- 
respond to two consecutive stochastic processes: A 
candidate R’ for (i + 1)-th sample is proposed accord- 
ing to the probability T(R’ + R;) and this move is 
then either accepted with the probability A(R’ + Ri) 
or rejected with the probability 1 — A(R’ + R;). If 
the move is accepted, the new member of the sequence 
is Risy = R’, otherwise it is Ri41 = R;. The length 
of the chain is thus incremented in either case. The 
acceptance probability A(R’ + R;), complementing 
some given T(R’ + Rı) and p(R) such that the sta- 
tionarity condition is fulfilled, is not unique. The 
choice corresponding to the Metropolis algorithm reads 


T(R & R’) eR’) 
' T(R — R) o(R) 


and depends only on ratios of T and p. Consequently, 
normalization of the trial wave function |W) is com- 
pletely irrelevant for the Monte Carlo evaluation of the 
quantum-mechanical expectation values. The freedom 
to choose the proposal probability T(R’ + R;) can 
be exploited to improve ergodicity of the sampling, for 
instance, to make it easier to overcome a barrier of 
low probability density separating two high-density 
regions. A generic choice for T(R’ + R;) is a Gaus- 
sian distribution centered at R; with its width tuned 
to optimize the efficiency of the sampling. 

The variational energy Eymc is a stochastic vari- 
able, and an appropriate characterization of the ran- 
dom error Eyyc — Eyz2, is thus an integral part of the 
variational Monte Carlo method. When the sampled 
local energies E (R;) are sufficiently well behaved [42], 
this error can be represented by the variance of Eymc. 
In such cases, the error scales as N~!/? and is pro- 
portional to fluctuations of the local energy. Reliable 


A(R’ + R) = min|1 


(5) 


estimation of the variance of Eymc is a non-trivial 
affair since the random samples {R;} generated by 
means of the Markov chain are correlated. These cor- 
relations are not known a priori and depend on the 
particular form of the transition matrix M that varies 
from case to case. Nevertheless, it is possible to esti- 
mate the variance without detailed knowledge of the 
correlation properties of the chain with the aid of the 
so-called blocking method [43]. 


The fluctuations of the local energy Fr, are reduced 
as the trial wave function |W) approaches an eigen- 
state of the hamiltonian, and Ey, becomes a constant 
when |W) is an eigenstate. In particular, it is crucial 
to remove as many singularities from Ey, as possible 
by a proper choice of the trial function. Section [4.1] il- 
lustrates how it is achieved in the case of the Coulomb 
potential that is singular at particle coincidences. 


The total energy is not the only quantity of inter- 
est and evaluation of other ground-state expectation 
values is often desired. The formalism sketched so far 
remains unchanged, only the local energy is replaced 
by a local quantity AL (R) = [Avy] (R)/Vr(R) cor- 
responding to a general operator A. An important 
difference between Ay, and the local energy is that 
fluctuations of Ay do not vanish when |W) is an 
eigenstate of H. These fluctuations can severely im- 
pact the efficiency of the Monte Carlo integration in 
(Wp|A| Wy) /(Wp|Wr), and the random error can decay 
even slower than M712 [42]. The trial wave function 
cannot be altered to suppress the fluctuations in this 
case, but a modified operator A’ can often be con- 
structed such that (¥r|Â'|Yr) = (W_|A|r) while 
the fluctuations of Ar, are substantially reduced 
A8}. 


2.2 Diffusion Monte Carlo 


The accuracy of the variational Monte Carlo method is 
limited by the quality of the trial wave function |W). 
This limitation can be overcome with the aid of the 
projector methods. In particular, the diffusion Monte 
Carlo method [49151] employs an imaginary time 
evolution 


Up(t)) =exp(-[A - Ext) Vr), (©) 


where the energy offset Ey is introduced to main- 
tain the wave-function norm at a fixed value. Formal 
properties of (6) can be elucidated by expanding the 
trial function |W) in terms of the hamiltonian eigen- 


states (2), which readily yields 


N¥o(0) = exp(—[Eo - Bre) [100 eolx) 


+ ew, Unley]. 


n=1 


The ground state |Wo) is indeed reached in the limit of 
large ¢ as long as the trial function was not orthogonal 
to |Wo) from the beginning. The requirement of a 
finite norm of |Wp(t)) translates into a formula 

Eo = Jim Er(t) (8) 
that can be used to obtain the ground-state energy. An 
alternative approach is to evaluate the matrix element 
Ep ur = (Wp (t)|H|Vr) /(Wp(t)|Vr) that asymptot- 
ically coincides with the ground-state energy, since 
(Wo|A |W) /(WolWr) = (Wolf |Vo) /(Wo|Wo). The in- 
tegration in Fy, can be performed stochastically 
in analogy with the VMC method, 


a J VŽ (R, t)Ur(R) [Hr] (R) BNR 


(Wp(t)|Yr)  YT(R) 
1 N 
x Epmc = NW - EL(Ri), (9) 


where the individual samples R; are now drawn 
from a probability distribution defined as p(R,t) = 
VER, t)Wr(R)/(Un(t)|Wr). 


2.2.1 Fixed-node/fixed-phase approximation 


The Monte Carlo integration indicated in (9) is pos- 
sible only if p(R,t) is real-valued and positive. Since 
the hamiltonians we usually deal with are symmetric 
with respect to time reversal, the eigenfunctions can 
be chosen real. Unfortunately, many-electron wave 
functions must necessarily have alternating sign to 
comply with the fermionic antisymmetry. In general, 
the initial guess |W) will have different plus and mi- 
nus sign domains (also referred to as nodal pockets 
or nodal cells) than the sought for ground-state wave 
function |Wo), which results in changing sign of p(R, t). 
In certain special cases, the correct sign structure of 
the ground state can be deduced from symmetry con- 
siderations [52H54], but in a general interacting system 
the exact position of the boundary between the pos- 
itive and negative domains (the so-called fermionic 
node) is unknown and is determined by the quantum 
many-body physics [55]. A number of exact properties 
of the fermionic nodes have been discovered [56H59], 
but a lot remains to be done in order to transform this 


knowledge into constructive algorithms for the trial 
wave functions. 

The problem with the variable sign of p(?, t) can 
be circumvented by complementing the projection (6) 
with the so-called fixed-node constraint [13], 


Up(R,t)Vr(R) >0 forall Randallt. (10) 


Doing so, limz+.|Wp(t)) only approximates |W), 
since the projection cannot entirely reach the ground 
state if the initial wave function |W) does not possess 
the exact nodes. The total energy calculated with 
this fixed-node method represents an upper-bound es- 
timate of the true ground-state energy because the 
projection (6) is restricted to a subspace of the whole 
Hilbert space when the constraint is implemented 
[60H62]. The fixed-node approximation has proved 
very fruitful in quantum chemistry [63] [64] as well as 
for investigation of the electronic structure of solids 
as testified by the applications reviewed in section [5] 

In calculations of extended systems and espe- 
cially metals, it is beneficial to allow for boundary 
conditions that break the time-reversal symmetry, 
since it facilitates reduction of finite-size effects (sec- 
tion [3.1). The eigenfunctions are then complex-valued 
and a generalization of the fixed-node approxima- 
tion is required. The constraint is replaced with 
Wp(t) = |Wp(t)| el” , where yr is the phase of the 
trial wave function Yr = |W] e'?t [65]. The phase pr 
is held constant during the DMC simulation to guar- 
antee that p(R,t) stays non-negative for all R and t. 
Additionally, a complex trial wave function |W) causes 
the local energy Er; to be complex as well. The ap- 
propriate modification of the estimate for the total en- 
ergy (9 coinciding with the asymptotic value of E(t) 
then reads 


1 N 
Epmc = WN > Re[Ei(Ri)| . (11) 


Analogous to the fixed-node approximation, the fixed- 
phase method provides a variational upper-bound es- 
timate of the true ground-state energy. Moreover, 
the fixed-phase approximation reduces to the fixed- 
node approximation when applied to real-valued wave 
functions. 


2.2.2 Sampling the probability distribution 


The unnormalized probability distribution that we 
wish to sample in the fixed-phase DMC method, 


F(R, t) = Up(R, t)¥r(R) = eNotes 5 
1 


referred to as the mixed distribution, fulfills an equa- 
tion of motion 


== F(R, t) = 
— SV2F(R,t) +V-[vp(R) F(R, t)] 


+ f(R,t) [Re[EL(R)] =F +t) Ert) (13) 


that is derived by differentiating (6) and with 
respect to time, combining the resulting expressions 
and rearranging the terms. The drift velocity vp in- 
troduced in is defined as vp = V ln|®¥r| and 
V denotes the 3N-dimensional gradient with respect 
to R. The equation of motion is valid in this form 
only as long as the kinetic energy is the sole non-local 
operator in the hamiltonian. Strategies for inclusion 
of non-local pseudopotentials will be discussed later 
in section The case of the fixed-node approxi- 
mation is virtually identical to (13), except that the 
local energy is real by itself. The following discussion 
therefore applies to both methods. 


The time evolution of the mixed distribution 
f(R,t) can be written in the form of a convolution 


F(R,t) = fore R’,t)f(R!,0) BNR’, (14) 


where f(R,0) = |Wr(R)|? and the Green’s function 
G(R « R’,t) = (R|G(t)|R’) is a solution of (13) 
with the initial condition G(R «+ R’,0) = d(R—-R’). 
Making use of the Trotter-Suzuki formula [66] [67], 
the Green’s function is approximated by a product of 
short-time expressions, 


G(t) = [Gesa(r) Gai (7) Garise(7)] + O(r), (15) 
$$$, ama 


N 


Gst (T) 


where 7 denotes t/M and the exact solution of 
is approached as this time step goes to zero. Con- 
sequently, the DMC simulations should be repeated 
for several sizes of the time step and an extrapola- 
tion of the results to r — 0 should be performed in 
the end. For simplicity, we show in only the 
simplest Trotter-Suzuki decomposition which has a 
time step error proportional to T. Commonly used are 
higher order approximations whose errors scale as T? 
or T?. The three new Green’s functions constituting 
the short-time approximation Gs can be explicitly 


written as 
Garitt(R + R’,T) (16) 
= [1—7V-up(R')]6[R— R' — vp (R’)r] 
+ O(r’), 


Gaia (R — R', T) (17) 
(R- wy 


(2rr)3N/2 aR | 2T 


Gesja(R — R’,7) (18) 
= exp |-7(Re [EL(R)] - Bx(t))| S[R-R'], 


and correspond to the three non-commuting oper- 
ators from the right-hand side of (13) in the or- 
der: drift (V - [vp(R)e]), diffusion (—V?/2e) and 
growth/decay (eè[Re[EL(R)] — (1 + tô:)Er(t)|). The 
drift and diffusion Green’s functions preserve the nor- 
malization of f(R,t) whereas the growth/decay pro- 
cess does not. 

The factorization of the exact Green’s function into 
the product of the short-time terms forms the basis 
of the stochastic process that represents the diffusion 
Monte Carlo algorithm. First, M samples {R;} are 
drawn from the distribution f(R,0) = |Wr(R)|? just 
like in the VMC method. Subsequently, this set of 
walkers evolves such that it samples the mixed distri- 
bution f(R,t) at any later time t. The probability 
distribution is updated from time t to t+ 7 by multi- 
plication with the short-time Green’s function, 


Reds / GalR Rr) f(Ri,t) PNR’, (19) 


which translates into the following procedure per- 
formed on each walker in the population: 


(i) a drift move ARarift = Up(R’)r is proposed 


(ii) a diffusion move AR gif = x is proposed, where 
x is a vector of Gaussian random numbers with 
variance T and zero mean 


(iii) the increment ARaritt + ARaig is accepted 
if it complies with the fixed-node condition 
Wr(R’)Ur(R’ + ARarite + ARaig) > 0, other- 
wise the walker stays at its original position; 
moves attempting to cross the node occur only 
due to inaccuracy of the approximate Green’s 
function (15), and they vanish in the limit 7 > 0; 
the moves ARaritt +A Raig are accepted without 
any constraint in the fixed-phase method 


(iv) the growth/decay Green’s function Gg q is ap- 
plied; several algorithms devised for this purpose 
are outlined in the following paragraph 


(v) at this moment, the time step is finished and the 
simulation continues with another cycle starting 


back at (i). 


After the projection period is completed, the algorithm 
samples the desired ground-state mixed distribution 
and the quantities needed for evaluation of various 
expectation values can be collected in step (vp. 

At this point we return to a more detailed 
discussion of several algorithmic representations of 
the growth/decay Green’s function Gq needed in 


step (iv). 


e The most straightforward way is to assign 
a weight w to each walker. These weights 
are set to 1 during the VMC initialization 
of the walker population and the applica- 
tion of Gg/q then amounts to a multiplication 
w(t +T) = w(t)W(R), where the weighting fac- 
tor is 


W(R) = exp [-7(Re [EL(R)] - Er(t))| . (20) 


Consequently, the formula for calculation of the 
total energy is modified to 


N -1 N 
Epmc = (Su) Dwi Re[ F(R] (21) 


i=1 


and the walkers remain distributed according 
to |Yr(R)|? as in the VMC method. This algo- 
rithm is referred to as the pure diffusion Monte 
Carlo method [68] [69]. It is known to be intrinsi- 
cally unstable at large projection times where the 
signal-to-noise ratio deteriorates [70], but it is 
still useful for small quantum-chemical systems 


[AF73]. 


e The standard DMC algorithm fixes the weights 
to w = 1 and instead allows for stochastically 
fluctuating size of the walker population by 
branching walkers in regions with large weight- 
ing factor W (R) and by removing them from 
areas with small W(R). The copies from high- 
probability regions are treated as independent 
samples in the subsequent time steps. The time 
dependence of the energy offset Er(t) provides a 
population control mechanism that prevents the 
population from exploding or collapsing entirely 
[50 [74]. The branching/elimination algorithm is 
much more efficient in large many-body systems 
than the pure DMC method, although it also 
eventually reaches the limits of its applicability 
for a very large number of particles . 


e An alternative to the fluctuating population are 
various flavours of the stochastic reconfigura- 
tion [I5] [70] [Z6H78]. These algorithms comple- 
ment each branched walker with high weighting 
factor W(R) with one eliminated walker with 
small W(R), and therefore the total number of 
walkers remains constant. This pairing intro- 
duces slight correlations into the walker popu- 
lation that are comparable to those caused by 
the population control feedback of the standard 
branching/elimination algorithm [75]. Keeping 
the population size fixed is advantageous for 
load balancing in parallel computational envi- 
ronments, since the number of walkers can be 
a multiple of the number of computer nodes 
(CPUs) at all times during the simulation. 


The branching/elimination process interacts in a subtle 
way with the fixed-node constraint. Since the walk- 
ers are not allowed to cross the node, the branched 
and parent walkers always stay in the same nodal cell. 
If some of these cells are more favoured (that is, if 
they have a lower local energy on average), the walker 
population accumulates there and eventually vanishes 
from the less favoured cells. Such uneven distribution 
of samples would introduce a bias to the simulation. 
Fortunately, it does not happen, since all nodal cells 
of the ground-state wave functions are connected by 
particle permutations and are therefore equivalent, see 
the tiling theorem in [56]. For general excited states 
this theorem does not hold and the unwanted depopu- 
lation of some nodal cells can indeed be observed. The 
problem is absent from the fixed-phase method, since 
it contains no restriction on the walker propagation. 


The branching/elimination algorithm is just one 
of the options of dealing with the weights along the 
stochastic paths. Another possibility was introduced 
by Baroni and Moroni as the so-called reptation al- 
gorithm [79], which recasts the sampling of both the 
path in the configuration space and the weight into 
a straightforward Monte Carlo process, avoiding thus 
some of the disadvantages of the DMC algorithm. The 
reptation method has its own sources of inefficiencies 
which can be, however, significantly alleviated [80]. 


This concludes our presentation of the stochastic 
techniques that are used to simulate the projection 
operator introduced in (6). We would like to bring 
to the reader’s attention that the algorithm outlined 
in this section is rather rudimentary and illustrates 
only the general ideas. A number of important perfor- 
mance improvements are usually employed in practical 
simulations, see for instance for further details. 


2.2.8 General expectation values 


So far, only the total energy was discussed in con- 
nection with the DMC method. An expression anal- 
ogous to (9) can be written with any operator A 
in place of the hamiltonian H. The acquired quan- 
tity Ay pw, = (Up|A|Vr)/(Wp| Ur), called the mixed 
estimate, differs from the pure expectation value 
(Up|A|Wp)/(Wp|Vp) unless A commutes with the 
hamiltonian. In general, the error is proportional to 
the difference between |Y p) and |W). The bias can 
be reduced to the next order using the following ex- 


trapolation [7] [50] 


(Wp|A|Wp) _ (Wp|A|Wr) (Wr AlWr) 
(Up|Wp) (UW p| Wr) (Ur|Vr) 
Up Ur 


T o( V(b Up) ) ie 


Alternative methods that allow for a direct evaluation 
of the pure expectation values have been developed, 
such as the forward (or future) walking [82], 
the reptation quantum Monte Carlo [79] [83] [84], or the 
Hellman—Feynman operator sampling [85] [86]. Due 
to their certain limitations, these techniques do not 
fully replace the extrapolation (22)—they are usable 
only for local operators and the former two become 
computationally inefficient in large systems. 

The discussion of the random errors from the end 
of section [2.1] applies also to the diffusion Monte Carlo 
method, except that the serial correlations among the 
data produced in the successive steps of the DMC sim- 
ulations are typically larger than the correlations in 
the VMC data. Therefore, longer DMC runs are neces- 
sary to achieve equivalent suppression of the stochastic 
uncertainties of the calculated expectation values. 


y (Ur| Ur) 


2.2.4 Spin degrees of freedom 


The DMC method as outlined above samples only the 
spatial part of the wave function, and the spin degrees 
of freedom remain fixed during the whole simulation. 
This simplification follows from the assumption of a 
spin-independent hamiltonian that implies freezing of 
spins during the DMC projection (6). This is indeed 
the current state of the DMC methodology as applied 
to electronic-structure problems: In order to arrive at 
the correct spin state, a number of spin-restricted cal- 
culations are performed and the variational principle 
is employed to select the best ground state candidate 
among them. 

Fixing spin variables is not possible for spin- 
dependent hamiltonians, such as for those containing 


spin-orbital interactions, since they lead to a non- 
trivial coupling of different spin configurations. In 
fact, spin-dependent quantum Monte Carlo methods 
were developed for studies of nuclear matter. A variant 
of the Green’s function Monte Carlo method 
treats the spin degrees of freedom directly in their full 
state space. Since the number of spin configurations 
grows exponentially with the number of particles, this 
approach is limited to relatively small systems. More 
favourable scaling with the systems size offers the aux- 
iliary field diffusion Monte Carlo method that samples 
the spin variables stochastically by means of auxiliary 
fields introduced via the Hubbard-Stratonovich trans- 
formation [89] [90]. Recently, a version of the auxiliary 
field DMC method was used to investigate properties 
of the two-dimensional electron gas in presence of the 
Rashba spin-orbital coupling [QI]. 


2.3 Pseudopotentials 


The computational cost of all-electron QMC calcula- 
tions grows very rapidly with the atomic number Z 
of the elements constituting the simulated system. 
Theoretical analysis [92] [93] as well as practical cal- 
culations [94] indicate that the cost scales as Z5-576-5, 
Most of the computer time spent in these simulations 
is used for sampling of large energy fluctuations in the 
core region, which have very little effect on typical 
properties of interest, such as interatomic bonding 
and low-energy excitations. For investigations of these 
quantities it is convenient to replace the core electrons 
with accurate pseudopotentials. A sizeable library of 
norm-conserving pseudopotentials targeted specifically 
to applications of the QMC methods has been built 
over the years [95}{100]. 

Pseudopotentials substitute the ionic Coulomb po- 
tential with a modified expression, 

Z s 
airs > Vir) +Ww (23) 

where V(r) is a local term behaving asymptotically as 
—(Z — Zeore)/r with Zeore being the number of elim- 
inated core electrons. The operator W is non-local 
in the sense that electrons with different angular mo- 
menta experience different radial potentials. Explicitly, 
the matrix elements of the potential W associated with 
I-th atom in the system are 


Imax 


l 
(RIWTR') =) SO (ulim) 


i l=0 m=-l 


X Wralrir)ôlrir — rir) m| Pir), (24) 


where |/m) are angular momentum eigenstates, riz is 
the distance of an electron from the J-th nucleus and 


fiz is the associated direction r;z/rir. Functions W7; 
vanish for distances riz larger than some cut-off ra- 
dius re, and the sum >; therefore runs only over 
electrons that are sufficiently close to the particular 
nucleus. 

Evaluation of pseudopotentials in the VMC method 
is straightforward, despite the fact that the local en- 
ergy Ex, itself involves integrals. As can be inferred 
from the form of the matrix elements (24), these are 
two-dimensional integrals over surfaces of spheres cen- 
tered at the nuclei. The integration can be imple- 
mented with the aid of the Gaussian quadrature rules 
that employ favourably sparse meshes {102}. 

The use of non-local pseudopotentials in the fixed- 
node DMC method is more involved, since the sam- 
pling algorithm outlined in section [2.2.2] explicitly as- 
sumed that all potentials were local. Non-local hamil- 
tonian terms can be formally incorporated by introduc- 
ing an extra member into the Trotter break-up (15), 
namely 


Gyioc(R — R', T) 


— Vr(R) -rů 1 
T TIR) (Rie |R’) 
= 6(R -R — 2 ®) Ri wWir’) + 0(72), (25) 


V7(R’) 


where W now combines the non-local contributions 
from all atoms in the system. This alone is not the 
desired solution, since the term involving the matrix 
element of W does not have a fixed sign and thus 
cannot be interpreted as a transition probability. 

To circumvent this difficulty the so-called localiza- 
tion approximation has been proposed. It amounts to 
a replacement of the non-local operator in the hamil- 


tonian with a local expression [93] {102} 


Ww > mR) = A) 


(26) 
Consequently, the contributions from W are directly 
incorporated into the growth/decay Green’s func- 
tion and no complications with alternating sign 
arise. Unfortunately, the DMC method with this ap- 
proximation does not necessarily provide an upper- 
bound estimate for the ground-state energy. The cal- 
culated total energy Epmc is above the lowest eigen- 
value of the localized hamiltonian, which does not 
guarantee that it is also higher than the ground-state 
energy of the original hamiltonian H. The errors in 
the total energy incurred by the localization approx- 
imation are quadratic in the difference between the 
trial function |W) and the exact ground-state wave 
function (102|. The trial wave functions are usually 
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accurate enough for the localization error to be practi- 
cally insignificant and nearly all applications listed in 
section |5| utilize this approximation. 

A method that preserves the upper-bound property 
of Epmc was proposed in the context of the DMC algo- 
rithm developed for lattice models [104]. The non-local 
operator W is split into two parts, W =Wt+W-, 
such that W+ contains those matrix elements, for 
which (R|W|R’/)Wr(R)Vr(R’) is positive, and W- 
contains the elements, for which the expression is neg- 
ative. Only the W* part is localized in order to obtain 
the approximate hamiltonian, 


z a WPR) 
W > W+ T(R) (27) 
One can explicitly show that the lowest eigenvalue of 
this partially localized hamiltonian is an upper bound 
to the lowest eigenvalue of the original fully non-local 
hamiltonian {104}. Recently, a stochastic represen- 
tation of the non-local Green’s function corre- 
sponding to W~ was implemented also into the DMC 
method for continuous space [105]. Apart from the 
recovered upper-bound property, the new algorithm 
reduces fluctuations of the local energy for certain 
types of pseudopotentials. On the other hand, the 
time step error is in general larger [105] (106), since the 
distinct treatment of the W+ and W~ parts of the 
pseudopotential essentially corresponds to a Trotter 
splitting of the growth/decay Green’s function 
into two pieces. Very recently, a more accurate Trotter 
break-up and other modifications improving efficiency 
of this method have been proposed for both continuous 
and lattice DMC formulations [107]. 

The localization approximation is directly applica- 
ble also to the fixed-phase DMC method. Adaptation 
of the non-local moves representing W~ to cases in- 
volving complex wave functions has not been reported 
yet, nevertheless, the modifications required should be 
only minor. 


3 From a finite supercell to the 
thermodynamic limit 


Quantum Monte Carlo methods introduced in the pre- 
ceding chapter can be straightforwardly applied to 
physical systems of a finite size, such as atoms and 
clusters of atoms. To allow investigation of bulk prop- 
erties of solids, the algorithms described so far have to 
be complemented with additional techniques that re- 
duce the essentially infinitely many degrees of freedom 
into a problem of manageable proportions. 


3.1 Twist-averaged boundary conditions 


In approximations that model electrons in solids as an 
ensemble of independent (quasi-)particles, it is possible 
to map the whole infinite crystal onto a finite volume 
so that the thermodynamic limit becomes directly ac- 
cessible. Hamiltonians of such models are invariant 
with respect to separate translations of electrons by 
any lattice vector R, that is, for each 7 we can write 


H(r1,r2,...,ri + R,...) 


= A(ri,r2,...,Ti,---)- (28) 


This invariance allows to diagonalize the hamiltonian 
only in the primitive cell of the lattice and then use 
the translations to expand the eigenstates from there 
into the whole crystal. Unfortunately, the explicit 
two-body interactions that we are set out to keep in 
the hamiltonian break the symmetry (28). The only 
translation left is a simultaneous displacement of all 
electrons by a lattice vector, which is not enough to 
reach the thermodynamic limit with a finite set of 
degrees of freedom. 

To proceed further we introduce artificial trans- 
lational symmetries with the aid of the so-called 
supercell approximation that is widely used within 
the independent-particle methods to investigate non- 
periodic structures such as lattice defects. We select 
a supercell having a volume Qs that contains several 
(preferably many) primitive cells. The whole crystal 
is then reconstructed via translations of this large 
cell by supercell lattice vectors {Rg}, which are a 
subset of all lattice vectors {R}. Simultaneously, we 
divide the electrons in the solid into groups containing 
N = PayQs particles, where pay is the average electron 
density in the crystal. This partitioning is used to 
construct a model hamiltonian, where electrons within 
each group interact, whereas there are no interactions 
between the groups, 


N 
=Y Sear) + Vee(R®)|. (29) 


The vector RO) = (ri, bibs, r) denotes coordinates 
of the electrons belonging to the I-th group. Note 
that these electrons are not confined to any particular 
region in the crystal. The supercell hamiltonian Hs 
consists of single-particle terms h, which encompass ki- 
netic energy as well as ionic and all external potentials, 


and of an electron-electron contribution 


; 1 
Vee(R) = >> ——_ 
a [ri = ryl 
j<i 
1 1 
+) D (30) 
TL? pees Pe Re] 


The first term in represents the explicit Coulomb 
interaction among electrons in the N-member group, 
and the second term mimics the interactions with the 
electrons outside the group. The physical meaning 
of the latter term is as follows: A set of images is 
associated with each electron j, and these virtual par- 
ticles are placed at positions r; + Rs so that they 
create a regular lattice. The combination of all images 
has the same average charge density as the original 
crystal and thus represents a reasonable environment 
for the selected N electrons. Each electron i then 
interacts with the arrays of charges associated with 
the other electrons in the group as well as with its 
own images. Only one half of the interaction energy 
with images is included in (30), the other half belongs 
to the rest of the system and is distributed among the 
other members of the sum in (29). The model hamilto- 
nian Hyode) approaches the original fully interacting 
hamiltonian as N increases and a larger fraction of 
the interactions has the exact form. 

Hamiltonians H. s and H, model possess the symmetry 
described by with lattice translations R replaced 
with Rs. Consequently, the eigenfunctions of Hg are 
many-particle Bloch waves 


N 
Wka(R) = UKa(R) exp (x ; yn) , (31) 


where œa is a many-body analogue of the band in- 
dex and K is the crystal momentum [108] [109]. The 
wave functions of the form (31) can be found in the 
same way as the single-particle Bloch waves within 
the independent-particle methods—as solutions to a 
problem of N particles confined to a simulation cell 
defined by vectors Lı, Lə and L} belonging to the 
set {Rs} and giving Qs = |L - (Lə x L3)|. The dy- 
namics of this finite N-particle system is governed by 
the hamiltonian Hg complemented with the so-called 
twisted boundary conditions |110 

Ti+ Lm,..., TN) 


= V kalfis -Tis TN 


Wka(ri, sa 
Je Em, (32) 
where i = 1,...,N and m = 1,2,3. The indistin- 


guishability of electrons implies that the phase factor 
is the same irrespective of which electron is moved, 
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Figure 1: Deviation of the twist-averaged to- 
tal energy (35) from the exact thermodynamic 
limit Bo, A = [Es(N) — Eoo]/N, for a three- 
dimensional gas of non-interacting electrons with 
dispersion relation €k = k?/2 at density p corre- 
sponding to rs = [3/(4mp)|‘/? = 1. The dashed 
line represents the average asymptotic decay of 
A that behaves as N71, 


and therefore only a single K vector common to all 
electron coordinates is allowed in and (32). Once 
the quantum-mechanical problem in the finite simula- 
tion cell is solved, wave functions for the whole crystal 
can be constructed. Since there are no interactions 
between individual N-particle groups, these wave func- 
tions have the form of an antisymmetrized product of 
the Bloch functions (31), namely 


WK {ar ({r}) = All] Ukor (R'”)] . (33) 
I=1 


The indicated antisymmetrization is straightforward 
as long as all Ky in the product are different, because 
each factor Y x,a, then comes from a disjoint part 
of the Fock space. The total energy corresponding 
to the wave function with the extra restriction 
Ky Æ Ky reads 


Co 


EK, ¢K yar} = So (Tko lfs|Y Ka) (34) 
I=1 


and the lowest energy is obtained by setting a; = 0, 
that is, by selecting the ground state for each of the 
different boundary conditions. Although unlikely, it 
is possible that the true ground state of Ĥ nodol falls 
outside the constraint Kz Æ Kr. In those cases, the 
expression with a; = 0 is an upper-bound esti- 
mate of the actual ground-state energy of Himodel with 
a bias diminishing as N increases. Taking into account 
the continuous character of the momentum K in the 
infinite crystal and the fact that all possible boundary 
conditions are exhausted by all K vectors within 
the first Brillouin zone, the ground-state energy per 
simulation cell can be written as 


` Qs 


I d°K (Wxo|Hs|Vxo). (35) 
LBZ. 


The total energy as well as expectation values of other 
periodic operators are calculated as an average over 
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N 


all twisted boundary conditions {110} [111]. In practice, 
the integral in is approximated by a discrete sum. 
The larger the simulation cell the smaller number of 
K points is needed to represent the integral, since the 
first Brillouin zone of the simulation cell shrinks and 
the K-dependence of the integrand gets weaker with 
increasing Ns. 

Formula is almost identical to the expression 
used in supercell calculations within the independent- 
particle theories, the only difference is that the number 
of electrons at each K point is fixed to N. This con- 
straint is benign in the case of insulators where the 
number of occupied bands is indeed constant across 
the Brillouin zone. In metals, however, the number 
of occupied bands fluctuates from one K point to an- 
other, and therefore the average does not coincide 
with the exact thermodynamic limit. Figure [I] pro- 
vides an illustration of the residual error. In principle, 
it is possible to remove this error with the aid of the 
grand-canonical description of the simulation cell [110], 
but this concept is not straightforward to apply since 
the supercell is no longer charge neutral. 


3.2 Ewald formula 


Our definition of the simulation-cell hamiltonian Hg 
in the preceding section was only formal and deserves 
a further commentary. It turns out that Voe is not 
absolutely convergent, and therefore it does not un- 
ambiguously specify the interaction energy. In partic- 
ular, the seemingly periodic form of the sums in 
does not by itself imply the desired periodicity of the 
hamiltonian. However, enforcing this periodicity as 
an additional constraint makes the definition unique 
and the resulting quantity is known as the Ewald en- 
ergy VE ) It can be shown that the requirement of 
periodicity is equivalent to a particular boundary con- 
dition imposed on the electrostatic potential at infinity 
[113]. The peculiar convergence properties of the 
sums in are a manifestation of the long-range 


character of the Coulomb potential—the boundary of 
the sample is never irrelevant, no matter how large 
its volume is. Consequently, careful considerations are 
required in order to perform the thermodynamic limit 
correctly. 

For the purposes of practical evaluation in QMC 
simulations, the Ewald energy is written as 


A - lim (vat) - =) , (36) 


where Vg(r; — rj) stands for an electrostatic potential 
at r; induced by the charge at rj together with its 
images located at r; + Rg. An explicit formula for 
the Ewald pair potential Vg reads [112] [14] 


Va(ri — rj) 


1 
= ——_—___ erfe(x|r; — 
2 ory Re iia 


An 1 G2 . 
te Dy aE ae Gs in) 
Gs40 


T 
Qgk?2 , 


Tj = Rs) 


(37) 


where Gs are vectors reciprocal to Rs, exp(iGs - 
Rs) = 1, and « is an arbitrary positive constant that 
does not alter the value of Vg. The omission of the 
Gs = 0 term in the reciprocal sum corresponds to 
the removal of the homogeneous component from the 
potential Vg. When evaluating the total energy of 
a charge neutral crystal, these “background” contri- 
butions exactly cancel among the Ewald energies of 
electron-electron, electron-ion and ion-ion interactions. 

The important feature of the Ewald formula is 
the decomposition of the slowly converging Coulomb 
sum into two rapidly converging parts, one in real 
space and the other in reciprocal space. The break- 
up is not unique (not only due to the arbitrariness 
of «) and can be further optimized for computational 


efficiency [115] [I6]. 


3.3 Extrapolation to the thermodynamic 
limit 

The total energy per electron ey = Eg/N evaluated 

according to the outlined recipe still depends on the 

size of the simulation cell. These residual finite-size 

effects can be removed by an extrapolation: Energy ey 

is calculated in simulation cells of several sizes and 


an appropriate function e#*(N) = €o + f(N) is subse- 
quently fitted through the acquired data. In the end, 
€> is the desired energy per electron in the thermody- 
namic limit, where the error term f(N) by definition 
vanishes, limy— oo f(N) = 0. Experience with a wide 
range of systems [10] [TI7H119] indicates that as long 
as the integral over the Brillouin zone in is well 
converged, the finite-size data are well approximated 
by a smooth function f(N) < 0 dominated by a 1/N 
contribution] The size extrapolation is therefore quite 
straightforward, although often computationally ex- 
pensive due to the relatively slow decay of the error 
term. 

The origin and behaviour of the finite-size effects 
is a subject of ongoing investigations with the aim 
to find means of accelerating the convergence of the 
total energy and other expectation values to the ther- 
modynamic limit. Furthermore, understanding the 
dependence of the finite-size errors on various parame- 
ters, such as particle density, can reduce the number of 
explicit size extrapolations needed to obtain quantities 
of interest. In calculations of equations of state (sec- 
tion 5.3), for instance, it is then sufficient to perform 
the extrapolation only at selected few, instead of all, 
electron densities [I19]. 

It turns out that, in the twist-averaged expecta- 
tion values (Hg) calculated in finite simulation cells, 
both the hamiltonian and the wave function contain 
biases that contribute to the 1/N asymptotics of the 
error term f(N). It was argued [I3] that the 
slow converging parts of the hamiltonian reside in the 
exchange-correlation energy 


^ 1 
Vio = (PP) | ole Velr—r'}ple') drd (38) 
QsxKXs 
defined as a difference between the expectation value 
of the Ewald energy VE )y and the Hartree term 
that describes the interaction of the charge densities 
p(r) = (p(r)). The Hartree energy is found to con- 
verge rather rapidly with the size of the simulation 
cell. In systems with cubic and higher symmetry, the 
leading contribution to f(N) can be written as 


_ Vxe _ Vxe 
Fxo(N) = -N — fim ay 
2 
So iy (39) 


This formula involves the static structure factor 


S(Gs) = 5, |[(A(Gs)A(-Gs)) - |(a(@s))["] , (40) 


N 


1The finite-size scaling depends on the dimensionality of the problem and the 1/N dependence corresponds to three-dimensional 


samples considered here. 
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where (Gs) is a Fourier component of the density 
operator. The exact small-momentum asymptotics of 
the structure factor in Coulomb systems is given by the 
random phase approximation (RPA) and reads 
S(Gs) ~ G2, which ensures that the limit in is 
well defined. In systems with lowered symmetry and 
for less accurate approximations, as is the Hartree— 
Fock theory where S(Gg) ~ Gg, the expression for 
fxc(N) must be appropriately modified [123]. 

The random phase approximation provides insight 
also into the finite-size effects induced by restricting 
the wave function into a finite simulation cell. Ac- 
cording to the RPA, the many-body wave function in 
Coulomb systems factorizes as 


WR) = War Re-Zu- a 


j<i 


where W, .. contains only short-range correlations and 
the function u(r) decays as 1/r at large distances. 
Such long-range behaviour is not consistent with the 
boundary conditions (32), and a truncation of this tail 
is therefore unavoidable. The corresponding finite-size 
bias is most pronounced in the expectation value of 


the kinetic energy T = (T’) and contributes an error 


term [121] 
fr(N) == - lim = 


1 

== Ae a G3 u(Gs), (42) 
where u(Gs) ~ 1/G2, is a Fourier component of u(r). 
Assuming that we are able to evaluate the expres- 
sions and (42), we can decompose f(N) into 
parts fN) = fxo(N) + fr(N) + f"(N), where f'(N) 
is substantially smaller than f(V) and the size extrap- 
olation is therefore better controlled. In the case of 
the homogeneous electron gas, the RPA provides ana- 
lytic expressions for the small momentum behaviour 
of the required quantities S(Gs) ~ G%/(2w,) and 
u(Gg) ~ 4r/(GRwp), where wp = 1/40 N/Qzg is the 
plasma frequency. Subsequently, the individual error 

terms simplify to 


It can be demonstrated that these two contributions 
completely recover the 1/N part of f(N), and that 
the residual term f’(N) scales as ~ N~4/9 [T23]. 
Application of the derived finite-size corrections to 
simulations of realistic solids is less straightforward 
since reliable analytic results are not available. The 
small momentum asymptotics of S(Gs) and u(Gs) 
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have to be examined numerically. Sufficiently large 
simulation cells are needed for this purpose, since the 
smallest nonzero reciprocal vector available in a su- 
percell with volume Qg is Gs ~ Oo, Utilization of 
the kinetic energy correction within DMC sim- 
ulations is further complicated by the fact that the 
function u(r) is not given as an expectation value of 
an operator, and thus it is not clear how it could be 
extracted from the sampled mixed distribution Yp Ur. 
One has to rely on the trial wave function alone to 
correctly reproduce the long-range tail (41), which can 
be a challenging task in simulation cells containing a 
large number of electrons. 


The above analysis employs exact analytic formu- 
lae or quantum Monte Carlo simulations themselves 
to find corrections to the finite-size biases. Alterna- 
tively, it is possible to adopt a more heuristic approach 
and estimate the finite-size effects within an approx- 
imative method. For instance, the error term f(N) 


can be rewritten as f(N) = eo = eer) + f”(N), 


where P% and ELDA) are total energies per particle 


provided by the local density approximation with two 
distinct exchange-correlation functionals, and f” (N) is 
anticipated to be considerably smaller than f(N) [4]. 
The exchange-correlation functional corresponding to 
a is constructed from properties of the homoge- 
neous electron gas in the thermodynamic limit (in 


other words, it is the standard LDA functional), the 


functional leading to ee is based on the homoge- 


neous electron gas confined to the same finite supercell 
as the quantum system under investigation. The latter 
functional is not universal and needs to be found for 
each simulation cell separately at the cost of auxiliary 
simulations of the homogeneous Coulomb gas. 


3.4 An alternative model for Coulomb 
interaction energy 


The expression for the electron-electron interaction 
energy Ve has two properties: (i) it is periodic and 
(ii) corresponds to an actual, albeit artificial, system of 
point charges. Although the latter property is concep- 
tually convenient, it is not really necessary, and any 
periodic potential that exhibits the correct behaviour 
in the limit of the infinitely large simulation cell is le- 
gitimate. Relaxation of the constraint (ii) in favour of 
faster convergence of the total energy per particle ey 
to its thermodynamic limit was explored in a series 
of papers [125], where the so-called model 
periodic Coulomb (MPC) interaction was proposed. 


(E) 


The replacement for the Ewald energy Vie? reads 
K 1 
(MPC) = ee 44 
Ve =e ae (44) 
Jct 
1 3 
+ 5 vel, —r) ae zho d?r 
i Qs 
- J mee) + | a(n) p(n’) drd2r’ 
E lr — r'|m 
QXsxQXs 


where |r—r’|m = ming, |r—r’+Rs| stands for the so- 
called minimum image distance. The operator vo PO) 
is constructed in such a way that the Hartree part of 
its expectation value is the same as in the Ewald 
method, whereas the slowly converging contribution 
to the exchange-correlation energy is removed. There- 
fore, the MPC interaction is essentially equivalent to 
the Ewald formula (36) complemented with the a pos- 
teriori correction (39). Instead of the structure factor, 
it is the one-particle density p that has to be evaluated 
as an extra quantity in this case (unless it is known 
exactly as in a homogeneous system). The explicit 
presence of the density p in the hamiltonian is incon- 
venient for the DMC method where the local energy 
is needed from the start of the simulation, that is, 
before the density data could be accumulated. The 
situation can be remedied by replacing the unknown 
density p with an approximation p4. For instance, the 
one-particle density provided by DFT is usually quite 
accurate. The error introduced by this substitution is 
proportional to (p — pa)? and further diminishes as 
the simulation cell increases. The Ewald and MPC 
energies per particle are therefore the same in the 
thermodynamic limit even if the approximate charge 


density is used [123] [126]. 


4 ‘Trial wave functions 


Accurate trial wave functions are essential for success- 
ful applications of the quantum Monte Carlo methods. 
The quality of the employed wave functions influences 
the statistical efficiency of the simulations as well as 
the accuracy of the achieved results. Equally impor- 
tant, especially for investigations of extended systems, 
is the possibility to quickly compute the wave func- 
tion value and its derivatives (Vr and V?W-r), since 
these quantities usually represent the most computa- 
tionally costly part of the whole simulation. Compact 
expressions are therefore strongly preferred. 

A significant part of the construction of the trial 
wave functions is optimization of the variational pa- 


rameters introduced into the functional form repre- 
senting Wy. It is a non-trivial task since the number of 
parameters is often large, Wr depends non-linearly on 
them, and the quantity to be minimized (Eymc or the 
variance of the local energy) is a fluctuating number. 
Several powerful methods addressing these problems 
have been developed during the years [[27}130] and 
even hundreds of parameters can be optimized with 
good efficiency nowadays. 


4.1 Elementary properties 


Since our aim is the electronic structure, and the elec- 
trons are subject to the Pauli exclusion principle, our 
trial wave functions have to be antisymmetric with re- 
spect to pair-electron exchanges. We assume collinear 
spins that are independent of electron positions, and 
therefore the full wave function Yr can be factorized 


as 
= Ny! Ny! 
tessa HE § (1° 
g N! 2 


x Yr(OR) |C -Th (45) 
Ny N 


where S = (01,...,0n) is a vector consisting of 
N = N+ + N, spin variables. The sum runs over 
all distinct states of N} up-oriented and N, down- 
oriented spins. In the case of N4 = 2 and N = 1 the 
spin states are |titels), | titst2) and |totsti), and 
the corresponding CR combinations are {r1, r2, T3}, 
{r1, r3, r2} and {r2, r3, r1}. The spatial-only part Yr 
is antisymmetric with respect to exchanges of par- 
allel electrons and its symmetry with respect to ex- 
changes of antiparallel electrons is unrestricted. The 
sum in with the appropriate sign factors (—1)° 
represents the residual antisymmetrization for the an- 
tiparallel spins. 

Both Yr and Wr are normalized to unity and 
identity (Wp|A|Wr) = (Y_|A|Wr) holds for any 
spin-independent operator A since the spin states 
|C{t ... fT ...L}) are mutually orthonormal. There- 
fore, it is usually sufficient to consider only the spatial 
part Wr of the full many-body wave function in appli- 
cations of the VMC and DMC methods, and we limit 
our discussion to Yr from now onf] 

Our goal is for the local energy HUy/Vey to be 
very close to a hamiltonian eigenvalue and fluctuating 
as little as possible. In systems with charged parti- 
cles interacting via the Coulomb potentials, it requires 
that the kinetic energy proportional to V? Yr contains 


In fact, the DMC algorithm is defined only for the spatial part Wr, consult sections and Decomposition then 
provides a hint how to calculate expectation values of spin-dependent operators from the sampled mixed distribution Y Yr. 
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singularities which cancel the 1/r divergencies of the 
potential. This cancellation is vital for controlling the 
statistical uncertainties of the Monte Carlo estimate 
of the total energy and takes place when Kato cusp 
conditions are fulfilled [132]. 

At electron-electron coincidences, these conditions 
can be formulated with the aid of projections of the 
trial wave function Vy onto spherical harmonics Yim 
centered at the coincidence point, 


l,m 
vit ) (riz, Tem R \ {ri, rj}) 


= f RYE OU). (46) 
Tij JAn 

In this definition, the following notation has been in- 
troduced: rij = |ri;| = |r:—r;]| is the electron-electron 
distance, Q;; is the spherical angle characterizing ori- 
entation of the vector rij, and Tem. = (ri + 17;)/2 is 
the position of the center of mass of the electron pair. 
The cusp conditions can then be written as 


for unlike spins and 


1 owam™ 4 
li T >= 48 
re vt”) Orij 4 ( ) 


for like spins. The expressions and differ 
because Wy is an odd function with respect to rj; in 
the latter case, which implies we 0) o, 

Analogous cusps occur in all-electron calculations 
when electrons approach nuclei. Unless the s-wave 


component yE ) vanishes (for a general discussion 
see [132]), it can be shown that 


1 agho 
D yo Orit 


=-Z), (49) 


where rj; is the electron-nucleus distance and Zy is 
the charge of the nucleus. 

A convenient functional form that meets the speci- 
fied criteria is a product of an antisymmetric part Va 
and a positive symmetric expression exp(—Ucorr) [133], 


Vr(R) = VA (R) exp[—Ucorr(R)]- (50) 


The Jastrow correlation factor exp(—Ucorr) incor- 
porates the electron-electron cusp conditions 
and (48), and Y4 ensures the fermionic character of 
the wave function. The electron-ion cusps can be 
included in either Va or in the corre- 
lation factor [136]. In simulations of extended systems, 
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the antisymmetric part obeys the twisted boundary 
conditions (section [3.1) and exp(—Ucorr) is periodic 
at the boundaries of the simulation cell. We discuss 
the individual parts of the trial wave function in 
some detail next. 


4.2 Jastrow factor 


The majority of applications fits into a framework set 
by the expression 


Ucorr(R) = 5 Xoi (r;) + > Ugi0; (ri; rj) , (51) 
i j<i 

where the functions x and u take a specific 
parametrized form and can depend also 
on spins of the involved electrons as indicated by the 
indices ø € {t,{}. The inclusion of the uncorrelated 
one-body terms y is important especially if the trial 
wave function is optimized with a fixed antisymmetric 
part Ya [139]. The two-body terms u are 
typically simplified to 


do tose; (ri: 77) + >) tee Tig) 
j<i j<i 
+ 5 Ueen(|Pig|, Parl, Irar), (52) 
j<i,I 

where Uee is an expression corresponding to a ho- 
mogeneous system and the electron-electron-nucleus 
term Ueen takes into account the differences between 
the two-body correlations in high-density regions near 
nuclei and in low density regions far from them. Spin 
indices were dropped to simplify the notation. The 
Ueen Contribution can usually be short ranged in the 
|r;;| parameter, whereas the simpler wee term is prefer- 
ably long ranged in order to approximate the RPA 
asymptotics (41) as closely as allowed by the given 
simulation cell [10] [136]. Of course, limited compu- 
tational resources can (and often do) enforce further 
simplifying compromises. In simulations of homoge- 
neous fermion fluids (electron gas, 3He), on the other 
hand, even higher order correlations were successfully 
included: three-particle [9] as well as four- 
particle [144]. 


4.3 Slater—Jastrow wave function 


The simplest antisymmetric form that can be used in 
place of Y4 in is a product of two Slater determi- 
nants, 


wR (R) = All (rf)... vh, wh, 
x ALT (rt) ek (rh, )] 
= det [Yf (r})] det [y4 ( rH], (53) 


where yf and y}, are single-particle orbitals that cor- 
respond to spin-up respectively spin-down electronic 
states and W?(r7) in the arguments of the determi- 
nants stands for a No x No matrix A®;. In quantum- 
chemical applications, a common strategy to improve 
upon the Slater wave function is to use a linear com- 
bination of several determinants, 


WER) = Y ca det ul nlr?) det [t mtrt) : 


(54) 
These multi-determinantal expansions are mostly im- 
practical for simulations of solids, since the number of 
determinants required to describe the wave function to 
some fixed accuracy increases exponentially with the 
system size. One case where multiple determinants are 
vital even in extended systems are fixed-node DMC 
calculations of excitation energies, since adherence 
to proper symmetries is essential for validity of the 
corresponding variational theorem and the 
trial wave functions displaying the correct symmetry 
are not always representable by a single determinant. 
In these instances, however, the expansions are 
short. 

In simulation cells subject to the twisted boundary 
conditions specified by a supercell crystal momen- 
tum K, the one-particle orbitals ~?, are Bloch waves 
satisfying 


Vikm(T + La) = VKmI(T) ets ’ (55) 


where a = 1,2,3 and m is a band index in the super- 
cell. Since the simulation cell contains several primitive 
cells, the crystal has a higher translational symmetry 
than implied by and the orbitals can be conve- 
niently re-labeled using m — (k,m’), where k and m’ 
are a momentum and a band index defined with re- 
spect to the primitive cell. The Bloch waves fulfill 
also 

Vrem (T + la) = Ykrem (r)i, (56) 
where la are lattice vectors defining the primitive cell. 
Assuming that the supercell is built as La = Nala 
with integers na, the momenta k compatible with 
fall onto a regular mesh 


Ny lL 3 (l2 x l3) 
i J2 ls x l J3 l x lə 


. ng l3- (l x =) on 


with indices ja running from 0 to ng — 1. This set of k 
points corresponds to the Monkhorst—Pack grid [146] 
shifted by a vector K. 

A number of strategies have been devised to de- 
termine the optimal one-particle orbitals for use in 


nə lə- (l3 x l) 


the Slater—Jastrow wave functions, which certainly 
differ from the Hartree-Fock orbitals that minimize 
the variational energy only when Uecorr = 0. Ideally, 
the orbitals are parametrized by an expansion in a sat- 
urated basis with the expansion coefficients varied to 
minimize the VMC or DMC total energy. The stochas- 
tic noise and the computational demands of the DMC 
method make the minimization of Epc extremely 
inefficient in practice. The VMC optimization of the 
orbitals was successfully performed in atoms and small 
molecules of the first-row atoms [130} [147] [748], but 
the method is still too demanding for applications to 
solids. 


To avoid the large number of variational param- 
eters needed to describe the single-particle orbitals, 
another family of methods has been proposed. The 
orbitals in the Slater—Jastrow wave function are found 
as solutions to self-consistent-field equations that rep- 
resent a generalization of the Hartree-Fock theory to 
the presence of the Jastrow correlation factor [139 149} 
[151]. These methods were tested in atoms as well as 
in solids within the VMC framework [152]. Un- 
fortunately, the wave functions derived in this way did 
not lead to lower DMC energies compared to wave 
functions with orbitals from the Hartree-Fock theory 
or from the local density approximation [149] [153]. It 
is unclear, whether the lack of observed improvements 
in the fermionic nodal surfaces stems from insufficient 
flexibility of the employed correlation factors or from 
the fact that only applications to weakly correlated 
systems were considered so far. 


An even simpler approach is to introduce a para- 
metric dependence into the self-consistent-field equa- 
tions without a direct relation to the actual Jastrow 
factor. An example are the Kohn—Sham equations cor- 
responding to some exchange-correlation functional, in 
which it is possible to identify a parameter (or several 
parameters) measuring the degree of correlations in 
the system and thus mimicking, to a certain degree, 
the effect of the Jastrow factor. Particular hybrid func- 
tionals with variable admixture of the exact-exchange 
component [37] were successfully employed for this 
purpose in conjunction with the DMC optimization, 
so that the variations of the fermionic nodal struc- 
ture could be directly quantified [154}{157]. Sizeable 
improvements of the DMC total energy associated 
with the replacement of the Hartree-Fock (or LDA) 
orbitals with the orbitals provided by the optimal hy- 
brid functional were observed in compounds containing 
3d elements. 


Evaluation of the Slater determinants dominates 
the computational demands of large-scale Monte Carlo 
calculations, and it is therefore very important to con- 
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sider its implementation carefully. Notably, schemes 
combining a localized basis set (atom-centered Gaus- 
sians or splines [158] [159]) with a transformation of 
the single-particle orbitals into localized Wannier-like 
functions can achieve nearly linear scaling of the com- 
putational effort with the system size when applied to 


insulators [T60}{162]. 


4.4 Antisymmetric forms with pair 
correlations 


Apart from the Pauli exclusion principle, the Slater de- 
terminant does not incorporate any correlations among 
the electrons, since it is just an antisymmetrized form 
of a completely factorized function, that is, of a prod- 
uct of one-body orbitals. A better account for cor- 
relations can be achieved by wave functions built as 
the appropriate antisymmetrization of a product of 
two-body orbitals. The resulting antisymmetric forms 
are called pfaffians and can generally be written as 


[591 [163] 


Np 
what(R) = A II a a | 
m=1 
N-—2Np 


x TT ver(rin) |, (58) 


where Np is the number of correlated pairs, Np < N/2. 
The two-body orbitals øf} that couple unlike-spin elec- 
trons (singlet pairs) are symmetric, whereas ¢/1 and 
o% (triplet pairs) are antisymmetric functions. The 
inclusion of the one-body orbitals yZ allows for odd N 
or for only partially paired electrons. The pfaffian 
wave functions can be viewed as compacted forms of 
particular multi-determinantal expansions (54). 

An important representative of the functional 
form is the Bardeen—Cooper-Schrieffer (BCS) 
wave function [164] projected onto a fixed number 
of particles, which is obtained from by con- 
sidering a singlet pairing in an unpolarized system 
(Ny = N, = N/2) with all two-body orbitals identi- 
cal. In that case, the antisymmetrization reduces to a 
determinant [165] 


WRS(R) = det [o"(r}, r})] , (59) 


where (rt, ry) is to be understood as a N/2 x N/2 
matrix A;;. In the quantum-chemical literature, this 
functional form is also known as the antisymmetrized 
geminal power. Note that the form of the BCS wave 
function does not by itself imply formation of Cooper 
pairs and their condensation, since the determinant 
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in reduces to the Slater wave function when 
the pair orbital is taken in the form 


N/2 


Poorer (Th TH) = dover elie}. (60) 
n=1 


The BCS—Jastrow wave functions were employed in 
investigations of ultra cold atomic gases (section [5.8) 
[166] [167] as well as in calculations of the electronic 
structure of atoms and molecules {169}. 

Trial wave functions with triplet pairing among 
particles were suggested in the context of liquid *He 
already two decades ago [165] [170]. It was realized 
only recently that even in these cases the exponentially 
large number of terms constituting the pfaffian can be 
rearranged in a way that facilitates its evaluation in a 
polynomial time, and therefore allows application of 
the pfaffian—Jastrow trial wave functions in conjunc- 
tion with the VMC and DMC methods [163] [171]. 


4.5 Backflow coordinates 


Another way to further increase the variational free- 
dom of the antisymmetric part of the trial wave func- 
tion is the backflow transformation Va(R) > Va(¥), 
where the new collective coordinates ¥ are functions 
of the original electron positions R. The designation 
“backflow” originates from an intuitive picture of the 
correlated motion of particles introduced by Feynman 
to describe excitations in quantum fluids [73]. 

In order to illustrate what is the origin of 
such coordinates, let us consider homogeneous 
interacting fermions in a periodic box with a 
trial wave function of the Slater-Jastrow type, 
U_(R) = det [exp(ik - ri)] exp[),<, y(riy)]- The Jas- 
trow factor is optimized so that its laplacian cancels 
out the interactions as much as possible within the 
variational freedom. Applying the kinetic energy op- 
erator to the Slater-Jastrow product results in local 
energy of the form 


[HVr|(R) _ 
m —Sasl®) 
— (v In|det [exp(ik - ri)] I) ; (v 5 (0)) , (61) 


i<j 


where we can qualitatively characterize Evar(R) 
as a mildly varying function close to a con- 
stant while the second term represents a non- 
constant “spurious” contribution, which appears as 
a scalar product of two fluxes. Consider now 
the following modification of the Slater—Jastrow 


form, W~(R) = det[exp(ik - x;)] exp| Mie; ylrij)l, 


where the single-particle coordinates are modified as 
xi = ri +); (rij) with (0) = 0. One can show that 
with a proper choice of the function ®(r;;), the lapla- 
cian of det |exp(ik-a;)] produces terms that cancel out 
most of the spurious contributions in the local energy 
given by (61). Of course, the backflow form generates 
also new non-constant terms so that the wave function 
has to be optimized for the overall maximum gain 
using variational strategies. 

In general, the new coordinates are written as 
x; = ri + €,(R) with € taken in a form analogous to 
the parametrization of the Jastrow factor Ucorr 
and (52). The vector € contains two-particle and pos- 
sibly higher order correlations and, in systems with ex- 
ternal potentials, also inhomogeneous one-body terms. 
The backflow transformation has been very success- 
ful in simulations and understanding of homogeneous 
quantum liquids [9] [144], and some progress 
has recently been reported in applying these techniques 
also to atoms and molecules [163] [174]. 


5 Applications 


In the last part of the article we go through selected 
applications of the quantum Monte Carlo methodology 
to the electronic structure of solids. In practically all 
listed cases, with the exception of sections [5.1] and 5-8] 
dealing with model calculations, the Slater—Jastrow 
functional form is employed as the trial wave function. 
The reviewed results therefore map out the accuracy 
that is achievable in the realistic solids when the mean- 
field topology of the fermionic nodes is assumed. It is 
shown that the quality of the DMC predictions is re- 
markable despite this relatively simple approximation. 


5.1 Properties of the homogeneous electron 
gas 


The homogeneous electron gas, also referred to as jel- 
lium, is one of the simplest many-body models that 
can describe certain properties of real solids, especially 
the alkali metals. At zero temperature, the model 
is characterized by the densities of spin-up and spin- 
down electrons, p} and p}, or, alternatively, by the 
total density p = p+ + py and the spin polarization 
Ç = |p — py|/p. It is convenient to express the den- 
sity p and other quantities in terms of a dimensionless 
parameter rs = [3/(47p)|!/3/ap, where ag is the Bohr 
radius. For example, the density of valence electrons 
in the sodium metal corresponds to rs % 4. 

The total energy of jellium is particularly simple 
since it includes only the kinetic energy of the electrons, 
the Coulomb electron-electron repulsion, and a con- 
stant which represents the interaction of the electrons 


with an inert uniformly distributed positive charge 
that maintains overall charge neutrality of the system. 
A straightforward dimensional analysis shows that the 
kinetic energy dominates the Coulomb interaction at 
high densities (small rs), where the electrons behave 
like a nearly ideal gas and the unpolarized state (¢ = 0) 
is the most stable. In the limit of very low densities, on 
the other hand, the kinetic energy becomes negligible 
and the electrons “freeze” into a Wigner crystal [189]. 


The homogeneous electron gas at zero temperature 
was one of the first applications of the variational and 
diffusion Monte Carlo methods. In the early inves- 
tigations [10] [I], only the unpolarized (¢ = 0) and 
fully polarized (¢ = 1) fluid phases were considered 
together with the Wigner crystal. Later, fluids with 
partial spin polarization were taken into account as 
well [90193]. The most accurate trial wave functions 
(the Slater—Jastrow form with backflow correlations) 
were used in reference [193] where it was found that 
the unpolarized fluid is stable below rs = 50 +2. At 
this density the gas undergoes a second-order phase 
transition into a partially polarized state, and the 
spin polarization ¢ then monotonically increases as 
the fluid is further diluted. Eventually, the Wigner 
crystallization density is reached, for which two DMC 
estimates exist: rs = 100 + 20 [II] and rs = 65 + 10 

192|. The discrepancy is presumably caused by the 

very small energy differences between the competing 
phases over a wide range of densities, and by uncer- 
tainties in the extrapolation to the thermodynamic 
limit. Advanced finite-size extrapolation methods, out- 
lined in section B]earlier, could possibly shed some new 
light on these quantitative differences. Indeed, recent 
calculations show further improvements in accuracy 
of the total and correlation energies [194]. A number 
of static properties of the liquid phases that provide a 
valuable insight into the details of the electron corre- 
lations in the jellium model and in Coulomb systems 
in general were evaluated by QMC methods as well 
[86] [TT] 193) 195]. 

The impact of the QMC calculations of the ho- 
mogeneous electron gas [IT] has been very significant 
because of the prominent position of the model as one 
of the simplest periodic many-body systems, and also 
due to the fact that the QMC correlation energy has 
become widely used as an input in density-functional 


calculations [35] [96]. 


The results quoted so far referred to the homo- 
geneous Coulomb gas in three dimensions. The two- 
dimensional gas, which is realized by confining elec- 
trons to a surface, interface or to a thin layer in a 
semiconductor heterostructure, has received similar 
if not even greater attention of QMC practitioners 
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Table 1: Cohesive energies of solids (in eV). Shown are 
DMC numbers unless only VMC data are available for 
the particular compound; those instances are marked 
with (*). The latest results are preferred in cases where 
multiple calculations exist. If not indicated otherwise, 
the experimental cohesive energies are deduced from the 
room-temperature formation enthalpies quoted in [188]. 


[197202]. The Wigner crystallization was pre- 
dicted to occur at rs = 37 +5 aF] a value that 
coincides with the density rs = 35 + 1 where a metal- 
insulator transition was experimentally observed later 
203]. 


5.2 Cohesive energies of solids 


The cohesive energy measures the strength of the chem- 
ical bonds holding the crystal together. It is defined as 
the difference between the energy of a dilute gas of the 
constituent atoms or molecules and the energy of the 
solid. Calculation of the cohesive energy is a stringent 
test of the theory, since it has to accurately describe 
two different systems with very dissimilar electronic 
structure. 

The first real solids whose cohesive energies were 
evaluated by a QMC method were carbon and silicon 
in the diamond crystal structure [204]. These 
early VMC estimates were later refined with the DMC 
method [180] [182] 206]. The most accurate re- 
sults to date are shown in table [I| where we have 
compiled the cohesive energies of a number of com- 
pounds investigated with the quantum Monte Carlo 
methods. Corresponding experimental data are shown 
for comparison. The electronic total energy calculated 
in QMC simulations is not the only contribution to 
the cohesive energy of a crystal, and the zero-point 
and thermal motion of the nuclei has to be accounted 
for as well, especially in compounds containing light 


compound QMC experiment 
Li 1.09 + 0.05 1.65 
1.57 +0.01 (x) 
Na 1.14 + 0.01 1.11 
1.0221 + 0.0003 [126 
Mg 1.51 +0.01 1.52 
Al 3.23 + 0.08 (x) 3.43 
MgH2 6.84 +0.01 6.83 
BN 12.85 + 0.09 (*) 12.9 [179 
C (diamond) 7.346+ 0.006 [180] 7.37 [181] 
Si 4.62+0.01 [182] 4.62 [183] 
Ge 3.85+0.10 {109} 3.86 
GaAs 4.9+0.2 (*) 6.7 
MnO 9.29 + 0.04 9.5 
FeO 9.66 + 0.04 9.7 
NiO 9.442+0.002 [186 9.5 
BaTiOs 31.2+0.3 187 31.57 


atoms. We refer the reader to the original references 
for details of these corrections. At present, a direct 
QMC determination of the phonon spectrum is gen- 
erally not practicable due to unresolved issues with 
reliable and efficient calculation of forces acting on the 
nuclei [48]. The effects due to the nuclear motion are 
thus typically estimated within the density-functional 
theory. 


Overall, the agreement of the DMC results with 
experiments is excellent; the errors are smaller than 
0.1 eV most of the time, including the Na and Mg 
elemental metals where coping with the finite-size ef- 
fects is more difficult. Notably, the diffusion Monte 
Carlo performs (almost) equally well in strongly cor- 
related solids represented in table [I] by 3d transition 
metal oxides MnO, FeO and NiO. The GaAs result is 
an obvious outlier with a systematic error of almost 
2 eV that the authors identify with the deficiencies 
of their pseudopotentials [184]. The two decades old 
application of the DMC method to the Li metal [175] 
is the only all-electron simulation in the list and its 
comparison to a subsequent pseudopotential calcula- 
tion [176] suggests that a large part of the discrepancy 
with the experiment is due to the fixed-node errors 
in the high-density core regions. It is likely that a 
substantial improvement would be observed if the all- 
electron calculations were revisited with today’s state 
of the art trial wave functions. 


3Note that in two dimensions the dimensionless density parameter rs is defined as rs = 1/(ap./7). 
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Table 2: Equilibrium lattice constants 


A3 
compound ao (A) Vo (A*) Bo (GPa) ao, equilibrium volumes Vo (per formula 
Li 3.556 + 0.005 1342 unit) and bulk moduli Bo for a number of 
3.477 [207] 12.8 [208] solids investigated with the QMC meth- 
ods. The first line for each compound 
Al 3.970 + 0.014 [177] 6517 contains QMC predictions, the second 
4.022 177 81.3 line shows experimental data. Theoret- 
GaAs 5.66 + 0.05 79+ 10 ical results for Li, Al and GaAs come 
5.642 77+1 from VMC simulations, the rest of the 
table corresponds to the DMC method. 
LiH 4.006 35.7 + 0.1 
4.07 +0.01 [212 32.2 + 0.03 
BN 11.812 + 0.008 [135] 378 + [135] 
11.812 + 0.001 395 + 2 213 
Mg 23.61+0.04 [106] 31.2+2.4 [106] 
23.24 214| 36.8+ 3.0 
MgO 4.23 216 158 
4.213 [217] 160 + 2 [217] 
MgH2 30.58 +0.06 (106| 39.5 1.7 
30.49 2185| — 
C 3.575 + 0.002 [219] 437 + [219] 
diamond) 3.567 442 + 4 
Si 5.439 + 0.003 [182] 103+10 {183} 
5.430 99.2 
SiO2 37.6 + 0.3 223| 32+6 
(quartz) 37.69 [224] 34 [224] 
FeO 4.324 + 0.006 [119 170 + 10 
4.334 225 = 180 226 


5.3 Equations of state 


The equilibrium volume Vo, the lattice constant ao, and 
the bulk modulus Bp = V(OE/OV)|v=v, are among 
the most basic parameters characterizing elastic prop- 
erties of a solid near the ambient conditions. Within 
QMC methods, these quantities are determined by eval- 
uating the total energy at several volumes around Vo 
and by fitting an appropriate model [227] [228] of the 
equation of state E(V) through the acquired data, 
see figure [2] for an illustration. Results of this proce- 
dure for a wide range of solids are shown in table [2] 
together with the corresponding experimental data. 
As in the case of the cohesion energy discussed in the 
preceding section, the raw QMC numbers correspond 
to the static lattice and corrections due to the motion 
of nuclei may be needed to facilitate the comparison 
with experimentally measured quantities. Particular 
details about applied adjustments can be inspected in 
the original papers. 

The data in table [2] demonstrate that the equilib- 
rium geometries predicted by the QMC simulations are 


very good and all lie within 2 % from the experiments, 
in many cases within only a few tenths of a percent. 
The agreement is slightly worse for the bulk moduli, 
where errors of several percent are common and in 
a few instances the mean values of the Monte Carlo 
estimates deviate from the experimental numbers by 
more than 10 %. Note, however, that determination of 
the curvature of L(V) near its minimum is impeded by 
the stochastic noise of the QMC energies and that the 
error bars on the less favourable results are relatively 
large. 

Quantum Monte Carlo methods are not limited 
to the covalent solids listed in table [2] Investigation 
of the equation of state of solid neon [229] represents 
an application to a crystal bound by van der Waals 
forces. Although the shallowness of the minimum 
of the energy—volume curve in combination with the 
Monte Carlo noise did not allow to determine the 
lattice constant and the bulk modulus to a sufficient 
accuracy, the DMC equation of state was still substan- 
tially better than results obtained with LDA and GGA. 
This example together with a recent study of interlayer 
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Figure 2: DMC total energies of the rock-salt 
(squares) and the NiAs (circles) phases of FeO. 
Statistical error bars are smaller than the symbol 
sizes. Lines are fits with the Murnaghan equa- 
tion of state. Inset: Difference between the Gibbs 
potentials of the two phases at T = 0 K; its in- 
tercept with the x axis determines the transition 
pressure P,. Adopted from [119]. 


E (eV/FeO) 


binding in graphite [230] illustrate that the diffusion 
Monte Carlo method provides a fair description of dis- 
persive forces already with the simple nodal structure 
defined by the single-determinantal Slater—Jastrow 
wave function. More accurate trial wave functions 
incorporating backflow correlations were employed to 
study van der Waals interactions between idealized 
metallic sheets and wires [231]. 

Calculations of the equations of state are by no 
means restricted to the vicinity of the equilibrium vol- 
ume, and many of the references quoted in table 
study the materials up to very high pressures. Such 
investigations are stimulated by open problems from 
Earth and planetary science as well as from other areas 
of materials physics. Combination of the equation of 
state with the pressure dependence of the Raman fre- 
quency {135} [219], both calculated from first principles 
with QMC, can provide a very accurate high-pressure 
calibration scale for use in experimental studies of 
condensed matter under extreme conditions [135]. 


5.4 Phase transitions 


Theoretical understanding of structural phase transi- 
tions often necessitates a highly accurate description 
of the involved crystalline phases. Simple approxi- 
mations are known to markedly fail in a number of 
instances due to significant changes in the bonding 
conditions across the transition. A classical example is 
the quartz-stishovite transition in silica (SiO2), where 
LDA performs very poorly and GGA is needed to 
reconcile the DFT picture with experimental findings 
232|. The diffusion Monte Carlo method has been 
employed to investigate pressure-induced phase tran- 
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sitions in Si [182], MgO [216], FeO and SiO2 
[223]. 

A transition from the diamond structure to the 
-tin phase in silicon was estimated to occur at 
P, = 16.5+0.5 GPa [182], which lies outside the range 
of experimentally determined values 10.3-12.5 GPa 
(see for compilation of experimental literature). 
Since the diamond structure is described very accu- 
rately with the DMC method as testified by the data 
in tables [I] and [2| it was suggested that the discrep- 
ancy is a manifestation of the fixed-node errors in the 
high-pressure 6-tin phase. This view is supported by 
a recent calculation utilizing the so-called phaseless 
auxiliary-field QMC, a projector Monte Carlo method 
that shows smaller biases related to the fermion sign 
problem in this particular case and predicts the tran- 
sition at 12.6 + 0.3 GPa [233]. It should be noted, 
however, that the volume at which the transition oc- 
curs was fixed to its experimental value in this later 
study, whereas the approach pursued in [182] was en- 
tirely parameter-free. 


In iron oxide (FeO), a transition from the rock- 
salt structure to a NiAs-based phase is experimentally 
observed to occur around 70 GPa at elevated tempera- 
tures [234] and to move to higher pressures exceeding 
100 GPa when the temperature is lowered [235]. The 
DMC simulations summarized in figure [2] place the 
transition at P, = 65+ 5 GPa [119]. This value repre- 
sents a significant improvement over LDA and GGA 
that both find the NiAs structure more stable than 
the rock-salt phase at all volumes. The agreement 
with experiments is nevertheless not entirely satisfac- 
tory, since the DMC prediction corresponds to zero 
temperature where experimental observations suggest 


stabilization of the rock-salt structure to higher pres- 
sures. Sizeable sensitivity of the transition pressure P, 
to the choice of the one-particle orbitals in the Slater— 
Jastrow trial wave function was demonstrated in a 
subsequent study [157], but those wave functions that 
provided higher P, also increased the total energies, 
and therefore represented poorer approximations of the 
electronic ground state. It remains to be determined, 
whether the discrepancy between the experiments and 
the DMC simulations is due to inaccuracies of the 
Slater—Jastrow nodes or if some physics not included 
in the investigation, for instance the inherently defec- 
tive nature of the real FeO crystals, plays a significant 
role. 

Investigations of phase transitions involving a liq- 
uid phase, such as melting, are considerably more 
involved due to a non-trivial motion of ions. An often 
pursued approach is a molecular dynamics simulation 
of ions subject to forces derived from the electronic 
ground state that is usually approximated within the 
density-functional theory. More accurate results would 
be achieved if the forces were calculated using quan- 
tum Monte Carlo methods instead. At present, this 
is generally not feasible due to excessive noise of the 
available force estimates |48|. Nevertheless, it was 
demonstrated that one can obtain an improved pic- 
ture of the energetics of the simulated system when 
its electronic energy is evaluated with the aid of a 
QMC method while still following the ion trajectories 
provided by the DFT forces 237]. 


5.5 Lattice defects 


The energetics of point defects substantially influences 
high-temperature properties of crystalline materials. 
Experimental investigations of the involved processes 
are relatively difficult, and it would be very helpful if 
the electronic structure theory could provide depend- 
able predictions. The role of electron correlations in 
point defects was investigated with the DMC method 
in silicon [206] [238] and in diamond [180]. The forma- 
tion energies of selected self-interstitials in silicon were 
found about 1.5 eV larger than in LDA, whereas the 
formation energy of vacancies in diamond came out 
as approximately 1 eV smaller than in LDA. These 
differences represent 20-30 % of the formation ener- 
gies and indicate that an improved account of electron 
correlations is necessary for accurate quantitative un- 
derstanding of these phenomena. 

Charged vacancies constituting the Schottky defect 
were investigated in MgO [239], and in this case the 
predictions of the DMC method differ only marginally 
from the results obtained within the local density 
approximation. The non-zero net charge of the su- 


percells employed in these simulations represents an 
additional technical challenge in the form of increased 
finite-size effects that require a modification of some of 
the size extrapolation techniques discussed in section [3] 


240) [241]. 


5.6 Surface phenomena 


Materials surfaces are fascinating systems from the 
point of view of electronic structure and correlation 
effects. The vacuum boundary condition provides 
surface atoms with more space to relax their posi- 
tions and surface electronic states enable the electronic 
structure to develop features which cannot form in 
the periodic bulk. This leads to a plethora of sur- 
face reconstruction possibilities with perhaps the most 
studied paradigmatic case of 7 x 7 Si(111) surface 
reconstruction. Seemingly, QMC methods should be 
straightforward to apply to these systems, similarly to 
the three-dimensional periodic solids. However, mainly 
technical reasons make such calculations quite difficult. 
There are basically two possibilities how to model a 
surface. One option is to use a two-dimensional pe- 
riodic slab geometry which requires certain minimal 
slab thickness in order to accurately represent the bulk 
environment for the surface layers on both sides. The 
resulting simulation cells end up quite large making 
many such simulations out of reach at present. The 
other option is to use a cluster with appropriate ter- 
mination that mimics the bonding pattern of the bulk 
atoms. This strategy assumes that the termination 
does not affect the surface geometries in a substantial 
manner. Moreover, it is applicable only to insulating 
systems. Given these difficulties, the QMC simulations 
of surfaces are rare and this research area awaits to 
be explored in future. 

The simplest possible model for investigation of 
surface physics is the surface of the homogeneous elec- 
tron gas that has been studied by DFT as well as QMC 
methods. The first QMC calculations [243] were later 
found to be biased due to complications arising from 
finite-size effects, especially due to different scaling of 
finite-size corrections for bulk and surface. Once these 
issues have been properly taken into account by Wood 
and coworkers [242], the QMC results have exhibited 
trends that were consistent with DFT and RPA meth- 
ods which are expected to perform reasonably well for 
this model system (see table 3). 

Applications to real materials surfaces are still very 
few. The cluster model was used in calculations of 
Si(001) surface by Healy et al. [244] with the goal 
of elucidating a long-standing puzzle in reconstruc- 
tion geometry of this surface, which exhibits regularly 
spaced rows of Si-Si dimers. The dimers could take 
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Table 3: Comparison of the surface energies (in erg cm~”) of the 
homogeneous electron gas calculated by a number of electronic struc- 
ture methods [242]. The DMC calculations were done with the LDA 2.07 


orbitals in the trial wave functions. 


two possible conformations: They can be either posi- 
tioned symmetrically or form an alternating buckling 
in a zig-zag fashion. While experiments suggested the 
buckled geometry as the low-temperature ground state, 
theoretical calculations produced conflicting results, 
in which various methods favoured one or the other 
structure. The QMC calculations [244] concluded that 
the buckled geometry is lower by about 0.2 eV/Sio. 
This problem was studied with QMC methodology 
also by Bokes and coworkers [245] who found that 
several systematic errors (such as uncertainty of ge- 
ometries in cluster models and pseudopotential biases) 
added to about 0.2 eV, and therefore prevented un- 
equivocal determination of the most stable geometry. 
This conclusion corroborated the experimental find- 
ings which suggested that at temperatures above 100 
K the distinct features of buckling were largely washed 
out and indicated that the effect is energetically very 
small. Very recently, the QMC study of this system 
has been repeated by Jordan and coworkers [246] with 
the conclusion that the buckled structure is lower by 
about 0.1 eV/Sig and that the highest level correlated 
basis set method which they used (CASPT3) is consis- 
tent with this finding. It was also clear that once the 
correlations were taken into account, the energy differ- 
ences between the competing surface reconstruction 
patterns were becoming very small. This brings the 
calculations closer to reality, where the two structures 
could be within a fraction of 0.1 eV/Sig as suggested 
by experiments. 


Perhaps the most realistic QMC calculations of 
surfaces have been done on LiH and MgO surfaces by 
the group of Alfé and Gillan who compared 
predictions of several DFT functionals with the fixed- 
node DMC method. The results showed significant 
differences between various exchange-correlation func- 
tionals. For the MgO(100) surface the best agreement 
with QMC results was found for the LDA functional, 
while for the LiH surface the closest agreement between 
QMC and DFT predictions was found for particular 
GGA functionals. 


Clearly, more applications are needed to assess the 
effectiveness of QMC approaches for investigation of 
surface physics. As we have already mentioned, the 
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Ts LDA GGA DMC RPA 

—608.2 —690.6 —563 +45 —517 
2.30 —104.0 —164.1 —82 + 27 —34 
2.66 170.6 133.0 179 + 13 216 
3.25 221.0 201.2 216 +8 248 
3.94 168.4 158.1 175 +8 182 


surfaces represent quite challenging systems for QMC 
methods. Nevertheless, we expect more applications 
to appear in the future since the field is very rich in va- 
riety of correlation effects that are difficult to capture 
by more conventional methods. 


5.7 Excited states 


The VMC and the fixed-node DMC methods both 
build on the variational principle, and they therefore 
seem to be applicable exclusively to the ground-state 
properties. Nevertheless, the variational principle can 
be symmetry constrained, in which case the algorithms 
search for the lowest lying eigenstate within the given 
symmetry class (provided, in the case of the DMC 
method, that the eigenstate is non-degenerate ), 
and thus enable access to selected excited states. 

Excitation energies in solids are calculated as differ- 
ences between the total energy obtained for the ground 
state and for the excited state. It is a computationally 
demanding procedure since the stochastic fluctuations 
of the total energies are proportional to the number of 
electrons in the simulation cell, whereas the excitation 
energy is an intensive quantity. Trial wave functions 
for excited states are formed by modifying the deter- 
minantal part of the ground-state Slater-Jastrow wave 
function such that an occupied orbital in the ground- 
state determinant is replaced by a virtual orbital. This 
substitution corresponds to an optical absorption ex- 
periment where an electron is excited from the valence 
band into the conduction band. The fact that both the 
original occupied orbital and the new virtual orbital 
necessarily belong to the same K point restricts the 
types of excitations that can be studied, since only a 
limited number of k points from the primitive cell fold 
to the given K point of the simulation cell, recall equa- 
tions (55)—(57). Clearly, the larger the simulation cell, 
the finer mapping out of excitations can be performed. 

Averaging over twisted boundary conditions (sec- 
tion [3.1) is not applicable to the calculations of the 
excitation energies, since both the ground state and 
the excited state are fixed to a single K point. This 
is not a significant issue, since finite-size effects tend 
to cancel very efficiently in the differences of the total 
energies calculated at the same K point. 


Table 4: Band gaps (in eV) of Mott insulators MnO and FeO cal- 


compound DMC experiment 
MnO 4.8 oa 0.2 3.9 a 0.4 
FeO 2.8+0.3 zx 2.4 [250 


DMC simulations following the outlined recipe were 
utilized to estimate the band gap in solid nitrogen 
and in transition-metal oxides FeO and 
MnO [248]. The gaps calculated for the two strongly 
correlated oxides are compared with experimental data 
in table [4| The ratio of the FeO and MnO gaps is 
reproduced quite well, but the DMC gaps themselves 
are somewhat overestimated, likely due to inaccuracies 
of the trial wave functions used for the excited states. 
A large number of excitations were calculated in sili- 
con [252] and in diamond [145] [253], and the obtained 
data were used to map, albeit sparsely, the entire 
band structure. In these weakly correlated solids the 
agreement with experiments is very good. Recently, a 
pressure-induced insulator—metal transition was inves- 
tigated in solid helium by calculating the evolution of 
the band gap with compression [254]. As illustrated in 
figure|3] (copyrighted material; not available in this ver- 
sion; see figure 1 in [254/), the DMC band gaps were 
found to practically coincide with the gaps calculated 
with the GW method. 


5.8 BCS—BEC crossover 


The repulsive Coulomb interaction considered so far 
is not the only source of non-trivial many-body ef- 
fects in the electronic structure of solids. A weak 
attractive interaction among electrons is responsible 
for a very fundamental phenomenon—the electronic 
states in the vicinity of the Fermi level rearrange into 
bosonic Cooper pairs that condense and give rise to 
superconductivity. The ground state of the system 
can be described by the BCS wave function Vgcg dis- 
cussed earlier in section [4.4] [164]. The Cooper pair 
is an entity that has a meaning only as a constituent 
of Vpcg. In order to form an isolated two-electron 
bound state, some minimal strength of the two-body 
potential is needed in three dimensions, whereas the 
Cooper instability itself occurs for arbitrarily weak 
attraction. When the interaction is very strong, the 
composite bosons formed as the two-electron bound 
states indeed exist and undergo the Bose-Einstein 
condensation (BEC). It turns out that a mean-field 
description of both the BCS and BEC limits leads to 
the same form of the many-body wave function, which 
indicates that the interacting fermionic system is likely 
to continuously evolve from one limit to the other when 
the interaction strength is gradually changed [256}{258]. 


culated with the fixed-node DMC method. Experimental data are 
provided for comparison. 


A large amount of research activity aimed at detailed 
understanding of this physics was stimulated by the 
possibility to realize the BCS-BEC crossover in exper- 
iments with optically trapped ultra-cold atoms [259]. 

In a dilute Fermi gas with short-ranged spherically 
symmetric inter-particle potentials, the interactions 
are fully characterized by a single parameter, the two- 
body scattering length a. The system is interpolated 
from the BCS regime to the BEC limit by varying 
1/a from —oo to oo. In experiments, this is achieved 
by tuning across the Feshbach resonance with the aid 
of an external magnetic field. Particularly intriguing 
is the quantum state of an unpolarized homogeneous 
gas at the resonance itself, where the scattering length 
diverges (1/a = 0). The only relevant length scale 
remaining in the problem in this case is the inverse 
of the Fermi wave vector 1/kp, and all ground-state 
properties should therefore be universal functions of 
the Fermi energy Ep. Since there is only a single 
length scale, the system is said to be in the unitary 
limit. The total energy can be written as 


3 
E = €Efree = E 5 ËF , (62) 


where Efree denotes the energy of a non-interacting 
system and € is a universal parameter. The universal- 
ity of € is illustrated in figure [4] that shows the ratio 
E/E'tvee as a function of the interaction strength cal- 
culated for three different particle densities using the 
diffusion Monte Carlo method with the trial wave func- 
tion of the BCS—Jastrow form. All three curves indeed 
intersect at 1/a = 0 with the parameter € estimated 
as 0.42 + 0.01 [166] (260}{262]. The energy calculated 
with the fermionic nodes fixed by the Slater—Jastrow 
wave function is considerably higher and would lead 
to € = 0.54 [166], which underlines the significance of 
particle pairing in this system. 

A further insight into the formation of the Cooper 
pairs is provided by evaluation of the condensate frac- 
tion that can be estimated from the off-diagonal long- 
range order occurring in the two-particle density ma- 
trix [263]. The condensate fraction a is given as 


N .. 
a=— lim par (r) (63) 


Too 


and the so-called projected two-particle density matrix 
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Figure 4: Fixed-node DMC energies of 38 
fermions in a cubic box with the periodic bound- 
ary conditions plotted as a function of interaction 
strength. BEC regime is on the left, BCS limit 
on the right. Shown are three particle densities p 
characterized by the dimensionless parameter rs 
defined in section [5.1] The simulations employed 
BCS—Jastrow trial wave function. Statistical er- 
ror bars are smaller than the symbol sizes. Data 
taken from |255]. 
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where rı corresponds to the spin-up state and r2 to 
the spin-down state. The evolution of a with interac- 
tion strength calculated with the DMC method [265] 
is shown in figure [5] (copyrighted material; not avail- 
able in this version; see figure 4 in [263)). It is found 
that approximately half of the particles participates in 
pairing in the unitary regime and this fraction quickly 
decreases towards the BCS limit, where only states in 
the immediate vicinity of the Fermi level contribute to 
the Cooper pair formation. Note that the condensate 
fraction vanishes if the Slater—Jastrow form is used in 
place of the trial wave function. 

The diffusion Monte Carlo simulations were used 
to study also the total energy and the particle density 
profile in the unitary Fermi gas subject to a harmonic 
confining potential [266}|268]. Due to the lowered 
symmetry compared to the homogeneous calculations 
referred above, the system sizes were more limited. To 
extrapolate the findings to a larger number of particles, 
a density functional theory fitted to the DMC data 
can be employed [269]. 


0.3 -0.2 -0.1 00 01 02 03 O04 05 
-1/(kp a) 
6 Concluding remarks 


In this article we have attempted to provide an 
overview of selected quantum Monte Carlo methods 
that facilitate calculation of various properties of cor- 
related quantum systems to a very high accuracy. Par- 
ticular attention has been paid to technical details 
pertaining to applications of the methodology to ex- 
tended systems such as bulk solids. We hope that we 
have been able to demonstrate that the QMC methods, 
thanks to their accuracy and a wide range of applica- 
bility, represent a powerful and valuable alternative to 
more traditional ab initio computational tools. 
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